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ADINC:  An  Implicit  Lagrangian 
Hydrodynamics  Code 

I.  INTRODUCTION 

This  paper  presents  and  describes  a new  software  package  called  ADINC  (ADiabatic  and 
INCompressible  flows)  which  is  designed  to  solve  the  mass,  momentum,  and  adiabatic  energy 
equations  for  a rather  general  one-dimensional  fluid  system.  To  permit  accurate  representation 
of  material  interfaces,  a fully  conservative,  Lagrangian  hydrodynamics  algorithm  has  been 
incorporated.  Implicit  time  differing  with  a quadratically  convergent  iteration  for  non-linear 
time  centering  permits  ADINC  to  take  long  timesteps  exceeding  the  Courant  sound  speed  limit. 
ADINC  is  designed  as  a user-oriented  package  with  a flexible  equation  of  state  which  may  vary 
from  cell  to  cell  and  provision  for  four  different  geometries.  This  write-up  not  only  describes 
the  techniques  employed  and  three  different  benchmark  test  problems,  it  serves  as  documenta- 
tion for  the  package.  To  these  ends,  listings  and  representative  outputs  are  included  in  the 
appendices. 

ADINC  was  written  to  circument  in  a single  software  package  a number  of  common 
numerical  dificulties  which  arise  in  the  simulation  of  many  different  fluid  dynamic  systems.  In 
some  fluids  the  physical  phenomena  under  study  vary  slowly  compared  to  the  time  it  takes 
sound  waves  to  cross  the  system.  Nevertheless,  substantial  compressions  and  expansions  occur 
so  the  incompressibility  assumption  is  invalid.  In  other  systems  rather  sharp  interfaces  between 
chemical  species  or  between  different  temperature  regions  must  be  maintained,  yet  the  fluids 
interact  dynamically  across  the  interface. 

The  first  class  of  problems  demands  an  implicit  treatment  of  sound  waves  in  general 
although  the  asymptotic  "slow-flow"  approach1  -J  works  well  in  many  cases.  In  many  combus- 
tion systems,  for  example,  heat  is  being  added  slowly  so  the  gas  expands  locally  where  the  heat 

is  being  released.  The  expansion  is  so  structured  that  it  guarantees  continued  spatial  constancy 
•Manuscript  submitted  April  20,  1979. 

1 


of  the  pressure.  Such  pressure  fluctuations  as  do  occur4  are  consistent  with  driving  the  fluid 
flows  which  bring  about  the  required  expansion.  In  these  slow  flow  systems  the  heating  and 
cooling  rates  from  chemical  reactions,  thermal  diffusion,  and  radiation  emission  and  absorption 
control  the  hydrodynamics  and  hence  determine  the  local  density  and  temperature.  The  pres- 
sure evolution  equation  can  be  solved  algebraically  for  the  divergence  of  the  velocity  at  any 
instant  of  time  while  the  vorticity  evolution  equation  is  integrated  to  advance  the  curl  of  the 
velocity.'  Given  the  curl,  the  divergence,  and  reasonable  boundary  conditions,  almost  every- 
thing we  need  to  know  about  the  flow  is  known.  When  the  pressure  fluctuations,  i.e.  the  sound 
waves,  themselves  are  required,  a stiff,  implicit  hyperbolic  wave  equation  results.  Such  equa- 
tions can  be  solved  numerically  but  in  multidimensions  the  cost  is  often  high  and  the  best  for- 
malism is  by  no  means  clear. 

In  the  study  of  ablation  and  deflagration  phenomena,  pressure  variations  are  important 
even  though  the  fluid  profiles  themselves  vary  slowly.  Here  again  the  sound  speed  in  the  high 
temperature  region  would  require  very  small  timesteps  if  an  explicit  integration  algorithm  is 
used.  Even  though  the  slow-flow  approximation  of  constant  pressure  is  not  valid  here,  the  gen- 
eralization to  implicit  hydrodynamics  is  valid-  as  long  as  shocks  are  absent  In  atmospheric 
flows,  in  turbulence  modelling  for  reactive  systems  and  in  studies  of  Rayleigh-Taylor  instabili- 
ties it  is  the  pressure  fluctuation  gradients  specifically  which  interact  with  the  local  density  gra- 
dients to  generate  vorticity.  Thus,  again,  some  form  of  implicit  hydrodynamics  or  implicit  pres- 
sure solution  is  required  in  these  and  similar  cases. 

In  imploding  liner  and  pellet  fluid  dynamics,  solid  and  higher  density  "fluids"  such  as 
compressed  I » • heavy  pusher  shells,  and  liquid  metal  rotating  liners  move  in  a way  which  is 
often  shock  free  and  which  generally  requires  non-irivial  equations  of  state.  In  LINUS6  8.  for 
example,  an  electrically  conducting  cylindrical  liquid  metal  shell  is  imploded  to  compress  a 
fusion  plasma  The  material  along  the  inner  surface  of  the  liner  starts  out  almost  perfectly 
imcompressible  but  compresses  by  as  much  as  50%  when  the  dynamic  pressure  from  geometric 
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convergence  reaches  its  maximum.  The  liquid  compression  is  important  and  hence  so  are  the 
detailed  dynamics  leading  up  to  that  compression.  The  sound  waves  in  the  liquid  metal  lead  to 
a "water  hammer"  effect7  8 which  may  have  severe  engineering  consequences. 

The  second  class  of  numerical  problems,  maintaining  sharp  interfaces  between  disparate 
materials  or  between  materials  in  widely  disparate  states,  often  demands  a Lagrangian  descrip- 
tion so  that  the  exact  interface  location  is  available  at  any  time.  Eulerian  computations  end  up 
. with  the  real  interface  spread  over  several  zones  from  the  numerical  diffusion  needed  for  stabil- 

ity. While  this  numerical  smearing  of  gradients  is  acceptable  in  many  cases,  in  many  others  it 
is  not  acceptable.  The  ADINC  package  uses  an  implicit,  fully  Lagrangian  algorithm  to  over- 
come these  two  classes  of  problems.  The  system  as  currently  modelled  is  assumed  to  have  no 
non-ideal  effects  such  as  viscosity.  Thus  energy  conserving  shocks  are  not  recovered  even 
though  large  amplitude  acoustic  waves  are  handled  accurately.  The  entire  package  is  structured 
so  that  rather  general  equations  of  state,  boundary  conditions,  and  flow  geometry  are  permitted 
Other  non-ideal  effects  can  be  introduced  to  the  calculation  via  time-step  splitting 

The  Courant  timestep  limit  is  overcome  by  using  an  implicit  finite  difference  algorithm 
with  adjustable  explicitness  coefficients  in  both  the  position  and  velocity  equations.  Since  the 
basic  equations  (Section  II)  are  nonlinear,  iteration  is  required  to  obtain  convergence.  ADINC, 
unlike  some  implicit  formulations,  continually  re-evaluates  the  nonlinear  terms  after  each  itera- 

ition  to  obtain  an  accurate  as  well  as  a stable  algorithm  when  changes  per  timestep  are  large 

I 

Since  this  increased  accuracy  has  the  computational  penalty  that  all  coefficients  for  the  iteration 
must  be  re-evaluated  each  iteration  cycle,  care  has  been  taken  to  develop  a quadratically  con- 
vergent iteration  encompassing  both  fluid  dynamics  and  equation  of  state  variations  simultane- 
ously. 

The  interface  resolution  problem  is  handled  by  making  the  calculation  fully  Lagrangian 
No  rezone  or  remap  capability  is  included  in  the  basic  package  so  numerical  diffusion  from  spa- 


tial differences  is  absent  from  the  model.  Future  additions  to  the  ADINC  package  will  include 
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a fully  Lagran^ian  adaptive  gridding  package  so  that  adequate  resolution  can  be  maintained  in 
regions  where  sharp  gradients  develop  or  where  zone  sizes  become  unacceptably  large  due  to 
fluid  expansion. 

Section  11  of  this  paper  discusses  the  basic  dynamical  equations  solved  by  the  ADINC 
package  and  the  various  choices  of  geometry  which  are  available.  Section  III  presents  the 
numerical  algorithms  used  to  integrate  these  equations  along  with  the  various  tricks  imple- 
mented to  maintain  accuracy  in  bizarre  situations.  The  timestep  conditions  built  into  the  pack- 
age are  also  discussed. 

Section  IV  is  devoted  to  the  structure  of  the  ADINC  package  and  reviews  the  three  major 
routines.  It  gives  calling  sequences  and  argument  lists  for  the  nine  entries.  ADINC  has  been 
constructed  as  a utility  package  to  advance  the  four  hydrodynamic  variables: 

r,  = position  (cm)  of  the  /-th  cell  interface 
V,  = velocity  (cm/sec)  of  the  /-th  all  interface 

p,  = density  (gm/cm3)  in  the  /-th  cell  between  interfaces  i- 1 and  /,  and 
P = pressure  (erg/cm3)  in  the  /-th  cell  between  interface  /-I  and  /'. 

ADINC  has  been  cast  into  a form  resembling  that  of  an  ordinary  differential  equations  package. 
The  user  requests  integration  of  the  system  of  equations  to  a certain  point  in  time  and  ADINC 
then  selects  the  number  and  length  of  timesteps  necessary  to  span  this  interval.  The  user  also 
has  control  of  various  error  and  integration  parameters  without  having  to  modify  the  ADINC 
code. 

Section  V tells  how  to  use  ADINC  but  a prospective  user  will  probably  find  Section  VI 
which  describes  three  test  problems  and  the  Appendices  equally  instructive.  The  three  test 
problems  were  designed  to  give  a prospective  user  a useful  background  of  experience  with  the 
program  in  different  regimes  and  with  different  problems.  The  Appendices  contain  listings  of 
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the  AD1NC  package  (Appendix  A)  and  a rather  general  lest  program  with  initialization  and  I/O 
routines  (Appendix  B).  The  present  version  of  ADINC  and  its  test  program  are  written 
entirely  in  64  bit  floating  point  arithmetic.  Thus  convergence  to  better  than  1 part  in  107  is  pos- 
sible for  problems  with  near  incompressibility  and/or  extreme  density  discontinuities  which 
require  this  accuracy.  A quadratically  convergent  algorithm  is  used  to  speed  joint  convergence 
of  the  nonlinear  fluid  dynamics  and  equation  of  state  physics.  The  only  portion  of  the  code  not 
appearing  in  the  appendices  are  the  vectorized  tridiagonal  solvers  and  these  are  documented  by 
Boris  in  NRL  Memorandum  Report  3408,  November,  1976. 

The  test  program  is  arranged  to  handle  up  to  five  distinct  layers  of  fluid  composed  of  up 
to  200  individual  finite  difference  cells.  Problems  in  one  of  four  geometries  can  be  set  up: 
Cartesian  coordinates,  cylindrical  coordinates,  spherical  coordinates,  and  a variable  power  series 
coordinates  for  treating  one-dimensional  nozzle-like  geometries.  The  boundary  conditions  are 
controlled  by  specifying  the  position  and  velocity  of  the  first  and  the  last  cell  interface.  Piston- 
like conditions  coupling  the  fluid  system  to  external  sources  and  sinks  of  energy  are  easy  to  set 
up. 

Appendices  C,  D,  and  E reprint  actual  output  from  the  three  test  problems  for  code 
verification  and  illustration  purposes.  ADINC  has  been  optimized  and  vectorized  for  the  Texas 
Instruments’  ASC  system  at  NRL. 


II.  THE  BASIC  EQUATIONS 


ADINC  solves  the  following  equations  for  mass  and  momentum  transport  in  one  dimen- 


sion: 
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I he  energy  evolution  equation  is  eliminated  by  using  an  adiabatic  equation  of  state  in  which  the 
entropy  is  assumed  constant  throughout  the  numerical  integration.  Non-adiabatic  processes 
such  as  external  heating,  thermal  conduction,  and  chemical  energy  release  can  be  added  to  Eqs. 
(1-2)  using  time  step  splitting  provided  sufficiently  short  timesteps  are  used  to  make  the  split- 
ting procedure  accurate.  In  the  reference  version  of  ADINC,  reproduced  as  Appendix  A of  this 
report,  the  equation  of  slate  of  the  fluid  in  each  cell  of  the  calculation  is 


p(p,  s.  ... ) -p,  . + (P/S)'!y.  (3) 

This  equation  of  state  with  p,  = 0 is  correct  for  adiabatic  compression  and  expansion  of  an 
ideal  gas.  In  that  case  1.2  < y < 1.7.  When  pc  ^ 0,  Eq.  (3)  gives  an  adequate  representation 
of  a mildly  compressible  liquid.  Water,  for  example,  has  p(.  = 1 gm/cc  and  S = 2.5  x 10"  in 
CGS  units.  Thus  in  this  crude  model  a pressure  of  250  Kbar  (2.5  x 10"  dynes/cm2)  causes  a 
substantial  increase  in  the  compression. 


During  an  ADINC  timestep  p,.,  y , and  5 are  treated  as  constants;  only  the  variation  of  p 
with  Pis  considered.  ADINC  does  not  use  the  temperature  T anywhere  and  uses  Eq.  (3)  in  the 
form  specified.  Rather  than  knowing  p and  asking  what  the  pressure  is,  ADINC  calculates  the 
fluid  density  given  an  approximation  to  the  pressure.  This  equation  of  state  density  is  com- 
pared to  the  density  derived  from  the  fluid  dynamics  via  Eq.  (1).  This  difference  is  iterated  to 
zero  using  a quadratically  convergent  implicit  solution  of  Eq.  (2)  which  delivers  an  improved 

pressure  approximation.  During  this  iteration  the  analytic  derivative  is  used  where  A is 
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the  volume  of  a computational  cell. 
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for  the  particular  equation  of  state  (3). 


(3) 


The  ADINC  package  is  written  in  a sufficiently  modular  form  that  replacement  of  Eq.  (3) 
with  another  equation  of  state  should  be  quite  straight  forward.  To  do  this  the  common  block 
/ADICOM/  would  have  to  be  modified  to  include  other  constants  of  the  fluid  motion  for  the 
various  materials  being  represented.  Thermochemistry  and  thermophysical  properties  of  realis- 
tic gases  could  be  included,  for  example,  so  the  effective  gas  constant  y can  be  made  to  display 
the  correct  variation  with  T.  As  another  example,  more  involved  equations  of  state  for  solid 
and  liquid  materials  can  be  included.  We  plan  to  use  Gardner’s  model9  in  studies  where  the 
transition  from  solid  to  plasma  must  be  treated  accurately. 


ADINC  uses  the  equation  of  state  in  the  form  p(P,  S,  . . .)  because  the  density  is  a much 
less  sensitive  function  of  the  pressure  than  the  pressure  is  of  the  density  for  liquids  and  solids. 
During  the  iteration  process,  finite  errors  in  pressure  and  density  are  expected.  In  the  other 
form  Pip,  S,  ...),  the  errors  in  density  p would  appear  as  wild  fluctuations  in  the  pressure. 
For  gases  and  plasma  the  two  forms  are  basically  of  ihe  same  accuracy.  There  is  a second 
related  reason  why  ADINC  uses  the  equation  of  state  in  the  case  p(P.  S,  ...),  a specific  form 
of  which  appears  in  Eq.  (3).  The  ADINC  package  is  specially  designed  to  deal  with  discontinui- 
ties in  zone  size  and  density.  When  a gas-solid  interface  is  encountered,  the  pressure  is  con- 
tinuous but  the  density  need  not  be.  Therefore  finite  differences  in  the  pressure  are  bound  to 
be  more  accurate  than  transfomed  differences  in  the  density. 


Equations  (1-3)  are  solved  in  the  form  shown  wi  tout  non-dimensionalization  or  scaling. 
The  package  is  designed  with  CGS  units  in  mind  but  these  appear  nowhere  explicitly  in  the  cal- 
culation. Therefore  non-dimensional  calculations  are  possible  without  modifying  the  program 
but  caveat  emptor. 
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To  complete  the  basic  equations  of  the  system,  the  geometry  of  the  calculation  must  be 
given.  ADINC,  as  presented,  allows  four  one-dimensional  systems  to  be  treated  by  selection  of 
the  integer  a Figure  1 shows  a general  one-dimensional  region  of  variable  cross  section  area 


Air).  The  volume  of  the  region  is  F(r)  = J A (r')  dr'.  The  numerical  algorithm,  discussed 

in  the  next  section,  uses  interface  areas  and  cell  volumes  exclusively  to  convey  geometry  infor- 
mation Therefore  substitution  of  other  ID  geometries  than  those  provided  for  should  be 
straightforward. 


a = 1,  Cartesian  coordinates 

Air)  = 1.0 
F(r)  = r 

a = 2,  Cylindrical  coordinates 

A (r ) = 2nr 
V(r)  = trr1 


(4.1) 


(4.2) 


a = 3,  Spherical  coordinates 

A ( r ) =*  Airr 2 
V(r)  =*  j nr3 

a = 4,  Power  series  (nozzle)  coordinates 

Air ) = Cr i -t-  G2r  + G}r 2 + G4r3  + G5r3 

.2  .4  5 

K(r)  = G,r  + ^5'5_ 


(4.3) 


(4.4) 


In  principle  the  volume  is  ambiguous  up  to  an  additive  constant,  the  volume  F(r|).  In 


practice  ADINC  deals  only  with  volume  differences  to  determine  the  incremental  volume  A of 


a computational  cell  so  this  constant  cancels  out  conveniently.  Using  the  symbol  A also  means 
that  confusion  with  the  velocity  Fcannot  occur. 


III.  THE  NUMERICAL  ALGORITHM 


Figure  1 shows  a schematic  diagram  of  the  computational  region  treated  by  the  ADINC 

package.  There  are  A' cells  of  volume  A,  (/  - 2,  3 N + 1)  bounded  by  /V  + 1 interfaces 

of  area  A,  (/  - I /V+  1).  The  interfaces  are  located  at  z,  (/  — 1 X + 1)  so 

A = A(r,)  where  A(r)  is  given  by  one  of  the  choices  a = 1,  2,  3.  4 in  Eq.  (4)  The  cell 
volumes  A,  = F(r,)  — V(r,_\)  are  the  difference  of  the  volume  integral  from  Eq.  (4)  at  the 
two  cell  interfaces.  In  the  development  to  follow  we  will  also  need  the  cell  center  locations 

R = (A,r,_ | + A,^,r,)/(A,  + A 

which  always  lies  between  and  r,  provided  the  interface  areas  are  positive. 

, In  ADINC  the  first  physical  cell  is  / — 2 and  it  lies  between  interfaces  r,  and  r2.  The  last 

physical  cell  is  / — X + 1 and  lies  between  interfaces  rv  and  rv+i-  Cells  1 and  N + 2 are  not 
used  currently  by  ADINC  but  have  been  left  available  for  future  use  in  complicated  boundary 
conditions  or  extrapolations.  Interface  1 is  the  left  hand  boundary  of  the  system  (r,  = rL)  and 
interface  N + 1 is  the  right  hand  boundary  of  the  system  (rv+i  = rR).  The  interfaces  are 
treated  in  a fully  Lagrangian  manner  and  therefore  the  interface  velocities  1 V\  are  defined  as 
well  as  the  interface  positions  {/-,}.  The  left  hand  boundary  velocity  (VL  = F,)  is  established 

I 

externally  as  is  the  right  hand  boundary  velocity  {VR  = F,).  The  interior  interface  locations 
and  velocities  are  the  quantities  which  ADINC  integrates  from  one  discrete  time  t to  the  next 
t + 8/  given  the  masses,  entropies,  and  other  cell  quantities  which  are  conserved  during  the 
motion. 


The  cell  interface  positions  (/•,)  satisfy 


dr, 

dt 


= V, 


(5) 


which  has  a straight  forward  discretization 


(5) 
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ADINC  GRID  STRUCTURE 
AND  VARIABLE  PLACEMENT 


Figure  I — Grid  structure  and  variable  definition  for  the  ADINC  package  Position  r and  Velocity  I are  defined  at 
cell  interfaces  while  density  p and  pressure  P are  defined  at  all  centers  Interface  area  A,  and  cell  volume  V are 
derived  from  the  instantaneous  interface  positions  (r,|. 


In  Eq.  (5  ) and  throughout  the  remainder  of  the  paper  the  superscript  "n"  indicates  variables  at 


the  "new"  time  t + 8 1 while  superscript  "o'  indicates  variables  at  the  "old"  time  t.  The  quantity 
t,  is  the  explicitness  parameter  for  the  interface  position,  0 ^ < 1.  When  «,  < 1,  the 
method  is  at  least  partially  implicit.  When  t — 1/2  the  method  is  centered  and  nominally  most 
accurate.  If  long  timesteps  are  comtemplated,  t < 1/2  is  required  for  Courant  stability  with 
strict  inequality  usually  required  to  deal  with  nonlinear  effects.  When  t,  — 0,  the  calculation  is 
fully  implicit,  i.e.  fully  forward  differenced.  This  case  is  most  stable  but  is  only  first  order  accu- 
rate. ADINC  uses  the  same  value  of  e,  at  every  interface  but  this  one  value  is  varied  from 
cycle  to  cycle.  Generalization  to  a spatially  varying  (*,)  is  quite  possible  but  left  for  the  future. 


The  momentum  equation  (2)  for  an  interface  velocity  is  as  follows: 


dK  s -1  dP 

dt  Pimerface* 


(6) 


where  the  density  and  pressure  gradient  are  needed  at  cell  interfaces.  The  discretization  used  in 
ADINC  is 


v°  - 


8r«, 


<p  8 r > 


8z ( 1 - Q 

<phr  >,  + |/2 


i + l/2 

(r+, 


(/?+.  - p°) 

- p?) 


(6) 


where  ev  is  the  explicitness  parameter  for  the  interface  velocity  and  has  the  same  properties 
described  above  for  tr.  The  quantities  «,  and  ev  are  distinct  in  ADINC  but  no  reason  has  been 
uncovered  to  date  for  using  different  values  in  an  actual  calculation.  The  interface  average 
indicated  as  <p8r>,+1/2  is  both  a spatial  and  temporal  average  as  described  below.  Physical 
considerations  are  used  to  define  <p8r>,  + )/2  so  the  discretization  in  Eq.  (6  ) is  insensitive  to 
numerical  errors  arising  from  large  density  discontinuities  at  the  interfaces. 


Figure  2 shows  two  cells  / and  / + 1 which  straddle  interface  i.  The  pressures  P,  and  P, 
are  defined  at  R,  and  /?/  + | as  shown  and  the  densities  p,  and  p,  + ( are  assumed  constant 
throughout  their  respective  cells.  Because  p,  and  p,  + , differ  spatially  (ignore  their  time 
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II 


ACCELERATION  MATCHING  ALGORITHM 


Figure  2 — Acceleration  matching  trick  for  density  discontinuities  at  interfaces  An  intermediate  interface  pressure  P' 
is  defined  such  that  the  acceleration  of  material  to  the  right  and  to  the  left  of  the  interface  is  matched  Most  the  pres- 
sure gradient  thus  appears  across  the  denser  fluid 

r 


i 

i 


variation  for  a moment),  the  straight  line  pressure  gradient  shown  would  impart  a different 
acceleration  to  the  fluid  just  to  the  right  and  to  the  left  of  interface  i.  If  the  fluid  were  permit- 
ted to  move  according  to  these  distinct  accelerations,  the  fluids  would  either  overlap  or  a gap 
would  open  up  at  interface  i after  a short  while.  To  prevent  this  a fictional  pressure  P’  is 
defined  at  interface  / such  that  the  acceleration  calculated  from  the  left  equals  the  acceleration 


calculated  from  the  right. 


p.+jr>\  + pj* 

+ n 


r = / 1 rn-  f<+'  = imr — 7T-  (8) 

PM,  - R ) p, +i(«,  + | - r,) 

In  terms  of  the  indicated  average  <p8r>,+1/2  in  Eq.  (6)  we  can  eliminate  P’  completely  from 
further  consideration  and  use 

<p8r >,+|/2  = p,  + |(/?,  + i - r ,)  + p,(r,  - R ,)  (9) 

to  define  the  spatial  part  of  the  free  average 


The  question  of  how  to  evaluate  the  average  (9)  in  time  arises  and  has  not  been  fully  set- 
tled The  major  points  to  consider  are  momentum  conservation,  nonlinear  instability  of  the 
overall  algorithm,  and  time-centering  accuracy.  Equation  (6),  when  multiplied  by  <p8r>,  + l/2 


and  summed,  yields 


£ K,"<p8r>  , + 1/2  - £ r,0<p8r  > ,+|/2  + boundary  terms 


where  the  telescoping  pressure  terms  cancel  except  at  the  boundary  of  the  computational 
region.  One  would  like  to  use  the  "old”  time  values  on  the  right  and  the  "new"  time  values  on 
the  left  to  give  a true  momentum  "integral"  — at  least  in  Cartesian  coordinates.  Since  the 
quantity  <p8r>,+|/2  is  conserved  in  Cartesian  Lagrangian  coordinates,  however,  it  doesn't 
really  matter  at  what  time  we  evaluate  <pf>r>,  + U2,  it  is  just  the  mass  associated  with  interface 
/.  In  non-Cartesian  coordinate  systems  the  momentum  integral  has  little  meaning  and  the 
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quantity  < pSr > , + 1/2  is  not  really  a constant  of  the  motion...  p, A,  is.  ADINC  uses  an  exactly 
time  centered  average  for  the  geometric  parts  of  <p8/->,  + 1/2.  The  density  which  appears  in 
the  expression  is  our  best  approximation  (latest  iteration)  to  the  new  density  for  reasons  of 
numerical  stability. 


Return  to  Eq.  (9).  ADINC  actually  uses 

<pbr  > , + \/2  = [p£,<**+,  - r,' 0 + p/V-  *,")]  (11) 

where  superscript  "P  stands  for  "previous"  and  indicates  the  latest  iterated  approximation  to  the 
"new”  value  of  the  variable,  in  this  case  The  superscript  "8"  is  used  to  indicate  the  exact 
"half’  time  average.  In  Eq.  (11) 

r*  = y (r,°  + rf) . Rj'  = y (/?/’  +/?,'’).  (12) 


We  make  no  representation  that  these  are  the  best  averages  or  that  extensive  testing  of  this 
aspect  of  ADINC  has  been  performed.  Perhaps  the  freedom  remaining  in  this  region  of  the 
calculation  can  be  used  to  further  improve  the  accuracy  and  veracity  of  the  algorithm.  We  do 
note  that  no  problems  arising  from  this  particular  choice  have  been  observed  to  date  in 
numerous  test  calculations  using  ADINC. 


Returning  to  the  complete  algorithm,  we  wish  to  find  a tridiagonal  equation  for  The 
momentum  equation  (6')  can  be  simplified  as  follows: 


where 


V?  = a,  - b,(P"+]  - P")  for  / = 2 N 


-^^7 iP:"  ~PX 

b,  = 8/(1  - «,.)/ <p8r>,  + |/j. 


The  equation  of  state  is  introduced  by  requiring  that  the  cell  volume,  computed  from 

the  equation  of  state  using  the  new  time  values  of  pressure,  equal  the  new  cell  volume  com- 
puted from  the  fluid  dynamics,  A,",rd).  At  any  iteration  Pthe  difference  is 


This  is  our  tridiagonal,  implicit,  linear  equation  for  the  estimated  new  pressures  [P"\  when  Eq. 
(13)  is  used  to  eliminate  1 K")  in  terms  of  (/?).  Expanding  this  out  for  completeness  gives 


Pieos) 

r~  + d;\a,  - b,(p,"+]  - />;)] 
-dr\a,- , - A, T-,)]  = c,. 


Equation  (21)  is  the  basic  equation  solved  by  AD1NC.  Iteration  is  necessary  because  Eqs. 
(17a')  and  (17b)  are  only  equalities  to  first  order.  The  iteration  is  quadratically  convergent, 
however,  because  only  second  order  terms  are  neglected  at  each  stage  of  the  iteration.  For  the 
interior  cells  / = 3,  4,  ....  N Eq.  (21)  may  be  written  as 

a;p^  + b;p:  + c;r+,  = d;  (22) 


A,  = -d,  b,^, 

. P(eos) 

R * j • 

" - Jp  , - c 

c;  = -d,+b,,  and 
A"  = c,  - d+a,  - d~ 


At  the  right  and  left  boundaries  of  the  system  the  new  velocities  are  given  externally  for  Eq. 
(19),  (19').  In  terms  of  the  tridiagonal  coefficients  of  Eq.  (23), 

D"i=D\ (above)  + d{V[ 

dn+\  =Av+t  (above)  - d„+ xVnR,  and  (24) 

A 2 = 0,  Cy+ 1 = 0. 

When  integrating  the  fluid  dynamic  equations,  ADINC  assumes  that  each  interface  moves 
in  a fully  Lagrangian  manner  according  to  Eq.  (5).  The  change  in  density  from  one  timestep  to 
the  next  in  a cell  is  therefore  given  simply  by  the  change  in  cell  volume  according  to  the  mass 
conservation  equation 

P"A,"  = \M,  = p,0A“  (25) 

When  individual  species  number  densities  must  be  followed,  they  are  also  advanced  by  Eq. 


(25). 


The  accuracy  and  stability  of  AD1NC  presupposes  the  monotonicity  of  the  interface  posi- 
tions {rj  with  increasing  / and.  of  course,  the  positivity  of  the  interface  areas  \A,).  Since  the 
algorithm  is  based  on  a discretization  of  the  continuum  fluid  dynamic  equations,  the  possibility 
exists  for  numerical  error  and  even  instability  caused  by  non-physical  crossing  of  cell  interfaces 
even  though  the  implicit  algorithm  given  above  is  nominally  stable  for  sound  waves  at  arbitrary 
timestep.  To  prevent  interface  crossings,  a Courant  condition  must  still  be  satisfied  for  the  flow 
velocities  1 V,}  even  though  | V\  « C5  = -J P/py  throughout  the  fluid  (gas).  In  the  reference 
versidn  of  AD1NC,  reproduced  in  Appendix  A,  the  maximum  timestep  which  still  prevents 
interfaces  crossing  adjacent  interface  positions  is  calculated  from  the  formula 


1 (rl+ , - r ,)  ( r , - r,_,) 

"al  2 Tv(  \V,\  ' \V\ 


where  a very  small  number  is  added  to  | I7, | to  prevent  dividing  by  zero. 


Equation  (26)  is  conservative  and  even  includes  the  factor  of  ~ in  case  two  interfaces  are 

moving  toward  each  other.  Generally  longer  timesteps  are  quite  acceptable.  By  rights  only  one 
of  the  two  terms  should  be  included  since  an  interface  can  be  moving  either  to  the  right  or  the 
left  but  not  both  at  the  same  time.  Furthermore,  interfaces  only  cross  due  to  a differential 
velocity,  not  an  absolute  one.  The  denominators  in  Eq.  (26)  should  be  | F,+1  — V\  and 
| V,  — F,_ 1 1 when  a net  motion  is  superimposed  on  relative  expansions  and  contractions. 


Situations  exist  where  the  timestep  calculated  from  Eq.  (26)  can  lead  to  trouble.  From 
Eq.  (5’)  one  can  see  that  the  cell  interface  positions  are  advanced  using  an  average  of  the  new 
and  old  velocities.  Since  the  timestep  has  to  be  estimated  at  the  beginning  of  a cycle,  it  is  quite 
possible  for  | F,°}  to  be  small  or  zero  while  | V")  can  be  large.  If  V”  becomes  large  enough,  r" 
can  cross  adjacent  interfaces  even  though  the  original  timestep  estimate  should  prevent  this. 
The  test  program  reproduced  in  Appendix  B also  contains  a user-supplied  maximum  timestep 


^max  10  provide  a limit  externally  because  no  fully  satisfactory  algorithm  has  yet  been 
developed  for  ADINC  to  estimate  the  limit  internally.  An  estimate  of  this  effect  would  require 
dealing  with  accelerations  as  well  as  velocities.  This  avenue  may  be  impractical  because  the 
acceleration  calculation  must  be  implicit  and  would  substantially  increase  the  computational  cost 
of  the  subroutine.  In  the  current  algorithm  non-convergence  of  the  pressure  iteration  by  a cer- 
tain amount  (discussed  below)  is  signalled  through  a diagnostic  message  printed  by  ADINC  but 
the  calculation  continues.  In  principle  non-convergence  after  a certain  number  of  iterations 
could  cause  re-initiation  of  the  computational  cycle  with  a smaller  timestep.  This  is  not 
included  in  the  current  version  for  three  reasons: 


1.  Extra  data  would  have  to  be  stored  to  re-initiate  the  calculation  at  additional  space 
and  time  expense. 

2.  The  ADINC  algorithm  already  subcycles  the  fluid  dynamic  calculation  when  the 
user  specifies  a longer  integration  interval  than  Eq.  (26)  permits.  Thus  the  logic 
is  already  quite  complicated. 

3.  Crossing  of  interfaces  is  usually  accompanied  by  rather  drastic  numerical  errors, 
not  nonconvergence,  therefore  recovery  is  problematical  at  best.  The  problem 
usually  blows  up  in  the  equation  of  state  calculation. 

For  these  reasons  I ve  chosen  to  put  the  onus  on  the  user  to  avoid  "sneaky"  interface  crossings. 
In  practice  the  timestep  given  by  Eq.  (26)  is,  in  fact,  conservative  not  overly  ambitious.  Few 
problems  have  arisen  to  date. 

ADINC  is  designed  to  be  used  as  a specific  partial  differential  equation  intergrator  in 
much  the  same  spirit  as  many  ordinary  differential  equation  packages  are  used.  It  attempts  to 
adjust  the  timestep  and  number  of  iterations  at  each  step  internally  to  maximize  speed  with 
good  accuracy.  Using  the  timestep  estimation  algorithm  given  above,  ADINC  subcycles  the 
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hydrodynamics  as  often  as  necessary  to  integrate  over  the  full  time  interval  which  the  user  has 
specified.  Since  the  opportunity  exists  for  an  inordinate  number  of  subcycles  if  something 
unexpected  happens  during  the  calculation,  the  number  of  subcycles  has  arbitrarily  been  limited 
to  100.  If  ADINC  performs  any  subcycles  at  all,  it  prints  a message  to  that  effect.  If  the  pack- 
age detects  that  more  than  100  subcycles  seem  to  be  required,  it  terminates  the  calculation 
printing  a suitable  diagnostic  message.  A user  may  wish  to  change  these  limits. 


Within  each  timestep  (or  subcycle)  the  package  performs  a convergence  calculation  on  the 
iterated  new  pressure  solution.  The  convergence  condition  used  for  each  cell  of  the  calculation 


< ,0-9  /4^ 


where  PmdK  is  the  largest  pressure  in  the  system,  Mm.M  is  the  largest  cell  mass  in  the  system, 
and  Mmm  is  the  smallest  cell  mass  in  the  system.  Equation  (27)  is  largely  heuristic  but  has 
been  found  to  work  extremely  well.  Using  the  timestep  condition  given  above,  convergence  is 
often  obtained  in  two  iterations  in  {/y}  and  almost  always  in  three  or  four.  1 he  maximum 
number  of  iterations  is  limited  to  six  because  this  should  be  sufficient  to  obtain  full  double  pre- 
cision accuracy  using  the  quadratically  convergent  algorithm  described. 


IV.  STRUCTURE  OF  THE  ADINC  PACKAGE 


I 

i 

j 


The  ADINC  package  consists  of  three  subroutines  containing  three  entries  each.  These 
FORTRAN  subroutines  are  vectorized  for  the  Texas  Instruments  ASC  system  at  the  Naval 
Research  Laboratory  and  are  reproduced  in  their  entirety  in  Appendix  A.  A rather  general  test 
program  and  several  simple  utility  subroutines  are  reproduced  as  Appendix  B.  Here  we  discuss 
the  three  major  routines  of  the  package  and  the  interactions  between  them.  The  next  section 
considers  more  practical  aspects  of  how  to  use  ADINC. 


The  three  routines  deal  with  the  geometry  of  the  problem,  the  equations  of  slate  of  the 
fluid,  and  the  fluid  dynamics  of  the  problem.  The  geometry  of  a problem  is  established  and 
controlled  by  the  subroutine  SETGEO  and  its  two  entries  USEGEO  and  DTFLOW.  The  calling 
sequences  and  arguments  for  these  entries  are  as  follows: 


CALL  SETGEO  (ALPHA.  GEOMCO) 
arguments: 

ALPHA  integer  I = Cartesian  coordinates 

2 = Cylindrical  coordinates 

3 = Spherical  coordinates 

4 = Power  Series  coordinates 

GEOMCO  real"8  array  contains  the  five  coefficients  G]  - G5of  Eq.  (4.4)  on  entry 
and  is  only  used  when  ALPHA  = 4 


CALL  USEGEO  (RAD 
arguments: 

RAD  real  *8  array 
AREA  real  *8  array 
RADC  real  ’8  array 

LAMC  real  '8  array 

N integer 


AREA.  RADC.  LAMC.  N) 

dimension  at  least  A+1  — contains  N+\  interface  positions  on  entry 
dimension  at  least  A+l  — contains  N+\  interface  areas  on  exit  (Eq. (4)  > 
dimension  at  least  N+2  — contains  ,V  celt  center  positions  on  exit 
in  locations  2,  . . . N+ 1 (Eq.  5)) 

dimension  at  least  A+2  — contains  A cell  volumes  on  exit 

in  locations  2.  . . V+  1 (Eq.  4)) 

contains  the  number  of  interior  cells  on  entry 


CALL  DTFLOW  IRAD.  VEL.  DTVAL.  N) 


arguments 
R AD  real  *8  array 
VEL  real  *8  array 
DTVAL  real  *8 
N integer 


dimension  at  least  .V+2  — contains  A+l  interface  positions  on  entry 
dimension  at  least  AM- 2 — contains  N+  I interface  velocities  on  entry 
contains  estimated  max  timestep  to  prevent  numerical  error  (on  exit) 
contains  the  number  of  interior  cells  in  entry 


SETGEO  Is  used  to  establish  the  geometry  of  the  problem  at  the  beginning  of  a calculation. 
Thereafter  USEGEO  is  called  to  fill  the  geometric  arrays  AREA,  RADC  and  LAMC  every  time 
a new  set  of  interface  positions  is  obtained  during  the  Lagrangian  calculation.  The  routine 
remembers  the  value  of  ALPHA  (and  the  coefficients  GEOMCO  where  appropriate)  until  SET- 
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GEO  is  called  to  set  up  a new  problem. 


DTFLOW  is  used  to  estimate  a maximum  timeslep  based  on  Eq.  (26)  for  use  in  the  fluid 
dynamic  calculation  performed  by  AD1NC.  The  reference  version  presented  in  Appendix  A 
calculates  a timestep  which  prevents  a new  interface  position  r,  + f'8rfrom  crossing  one  of  the 
original  adjacent  interface  locations,  r, or  r,_ , . 


The  equation  of  state  of  the  fluids  in  each  cell  of  the  computation  is  calculated  by  the  sub- 
routine SETMAT  and  its  two  entries  SETEOS  and  USEEOS.  The  calling  sequences  and  argu- 
ments for  these  entries  are  as  follows: 


I 

! 


CALL  SETMAT  (MATER. M ASS.G  AMM  A.RHOCON.N) 
arguments 

MATER  real  '8  array  dimension  at  least  ,V+2  — contains  A/ cell  material  identifiers  (i-2.  . Sri1 

MASS  real  "8  array  dimension  at  least  ,V+2  — contains  cell  masses  on  entry 

GAMMA  real  '8  array  dimension  at  least  iV+2  — contains  on  entry  .V cell  y values  from  Eq.  (31 

RHOCON  real  *8  array  dimension  at  least  ,V+2  — contains  on  entry  N cell  p values  from  EQ  (3) 

N integer  contains  the  number  of  interior  cells  on  entry 

CALL.  SETEOS  (RHO.PRE.N) 
arguments. 

RHO  real  *8  array  dimension  at  least  ,V+2  contains  N cell  densities  on  entry 

PRE  real  ’8  array  dimension  ai  least  N+ 2 — contains  N pressures  on  entry 

,V  integer  contains  the  number  of  cells  on  entry 

CALL  USEEOS  (RHO.Pkfc  LAMEOS.DLAMDP.N) 
arguments 

RHO  real  *8  array  dimension  at  least  N+2  — contains  IV cell  densities  p(P,  S,  ) on  exit 

PRE  real  *8  array  dimension  at  least  N+2  — contains  N pressures  on  entry 

LAMEOS  real  *8  array  dimensional  least  N+2  — containts  /V cell  volumes  on  exit 
01.  AMDP  real  *8  array  dimension  at  least  N+2  — containts  d \/dP  for  each  cell  on  exit 
N integer  contains  the  number  of  interior  cells  on  entry 


SETMAT  is  used  lo  place  the  equation  of  state  cell  constants  into  the  common  block  /ADI- 
COM/  which  communicates  equation  of  state  information  among  the  three  routines  (nine 
entries)  of  the  ADINC  package  These  constants  are  those  conserved  quantities  and  local  pro- 
perties of  the  fluid  cells  used  to  connect  the  grid  configuration  at  one  time  to  that  at  another 
adiabatically.  Entry  SETEOS  is  called  to  determine  the  entropy  for  given  values  of  cell  density, 
ceil  pressure  and  interface  positions.  Since  thermal  conduction,  external  healing,  chemical 
energy  release,  etc.  all  change  the  entrony  but  not  the  other  constants  in  the  equation  of  state 
a separate  entry  from  SETMAT  is  provided  to  reduce  the  cost  of  resetting  the  equation  of  state 
to  account  for  these  non-ideal  phenomena. 
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USEEOS  is  the  entry  provided  to  determine  the  cell  density  expected,  given  the  cell 
pressures  and  all  the  equation  of  state  constants  in  ADICOM,  Also  returned  to  the  user  is  the 
cell  volume  expected  and  the  rate  of  change  of  change  of  cell  volume  with  pressure.  These 
quantities  depend  on  the  equation  of  state  only  and  are  used  in  the  ADINC  iteration  which 
reconciles  the  cell  volumes  computed  from  the  fluid  dynamic  equations  with  those  needed  to 
satisfy  the  equation  of  state.  While  SETMAT  and  SETEOS  need  only  be  called  only  once  by 
the  user  in  many  types  of  calculations,  USEEOS  is  called  at  least  two  times  periteration  by 
ADINC  and  therefore  entails  the  lion's  share  of  the  computational  expense. 

The  common  block  /ADICOM/  is  declared  as  follows: 

PARAMETER  NPT  = 202 

INTERGER  NCELLS 

REAL  *8  MATERC(NPT).  MASSC(NPT).  GAMMAt(NPT) 

REAL  *8  ENTC(NPT).  RHOC(NPT) 

COMMON  /ADICOM/  MATERC,  MASSC.  GAMMAC.  ENTC,  RHOC.  NCELLS 

and  appears  in  the  three  ADINC  routines  in  exactly  the  same  form.  If  a user  wishes,  more 
complicated  formulae  may  be  used  for  an  equation  of  state.  The  additional  arrays  needed 
would  then  be  added  to  each  realization  of  /ADICOM/  and  the  appropriate  computations  to 
SETMAT,  SETEOS,  and  USEEOS.  The  version  provided  allows  a relatively  wide  selection  of 
problems  to  be  tackled  with  a relatively  simple  formulation  that  is  inexpensive  to  compute. 

The  quantity  MATERC(I)  is  the  material  and  shell  identification  of  cell  i.  It  is  a floating 
point  number  of  the  form  L.MM  where  0 < L < 10  is  the  layer  identifier  and  0 < MM  < 100 
is  the  local  material  identifier.  The  layer  number  is  used  by  the  diagnostics  to  perform  various 
partial  sums.  The  layer  can  contain  different  materials  from  cell  to  cell  but  the  layer  number 
must  increase  monotonically  from  cell  2 to  cell  A'+l.  The  material  identifier  MM  is  not  used 
in  the  current  version  but  is  available  for  flagging  which  of  several  computationally  distinct 
equations  of  state  are  to  be  used  in  the  future.  MASSC(l)  is  the  mass  (nominally  in  grams)  of 
cell  i.  GAMMAC(I),  RHOC(I),  and  ENTC(I)  are  the  equation  of  state  conserved  constants 
for  cell  / which  appear  in  Eq.  (3). 


Double  precision  is  used  throughout  the  floating  point  computations  so  that  problems  with 


large  disparity  in  cell  mass  or  with  nearly  incompressible  fluids  can  be  treated  aceuratel>  The 
basic  convergence  criterion  of  10  9 for  the  quadratically  convergent  iteration  in  ADINC  is  far 
more  accurate  than  single  precision  on  the  ASC  permits.  We  are  using  a 64  word  length  here 
with  an  8 bit  (hexidecimal)  exponent.  The  versions  of  the  equation  of  state  and  geometry  rou- 
tines are  also  fully  vectorized  for  efficient  computation  on  the  ASC. 


The  basic  fluid  dynamic  calculation  is  performed  by  ADINC  and  its  two  entries  AD1NCO 
and  SETEPS. 


CALL  ADINC  (RAD.VEL.RHO.PRE.N.DTIN.CYCLE.RRNEW.RLNEW.VRNEW.VLNEWI 


argumems: 

RAD  real  *8  array 
VEL  real  *8  array 
RHO  real  "8  array 
PRE  real  ‘8  array 


dimension  ai  least  N+  1 — contains  N+  1 interface  positions 
dimension  at  least  Ai+I  — contains  ,V+ I interface  velocities 
dimension  at  least  V+2  — contanis  .V  cell  densities 
dimension  at  least  /V+2  — contains  jV cel!  pressures 


These  four  arrays  of  physical  variables  contain  the  old  timestep  values  on  entry  and  will  be  filled  with  the 
new  timestep  values  DTIN  later  on  exit.  Users  should  store  old  values  in  the  external  program  il  desired 
The  test  program  shows  how  to  do  this 


N integer 
DTIN  real  *8 
CYCLE  integer 
RRNEW  real  *8 
RLNEW  real  '8 
VRNEW  real  *8 
VLNEW  real  *8 


contains  the  number  of  interior  cells  on  entry 
the  integration  timestep  on  entry 
the  timestep  number  for  identification 
on  entry  the  new  right  boundary  position 
on  entry  the  new  left  boundary  position 
on  entry  the  new  right  boundary  velocity 
on  entry  the  new  left  boundary  velocity 


The  boundary  interface  values  are  specified  at  the  end  of  the  desired  timestep  since  they  are  assumed  to  be 
driven  externally.  Usually  the  positions  are  fixed  and  the  velocties  zero  RAD(,Y+I),  RAD(l).  V'EL(.V+1). 
(VEL(I)  should  contain  the  corresponding  old  time  values. 


CALL  ADINCO  (MODE.  CYCLE) 
arguments 

MODE  integer  = 0 if  reinitializing  NCALL 

= I if  NCALL  cumulative-as  described  in  the  next  section 
CYCLE  integer  the  timestep  number  for  identification 


CALL  SETEPS  (ERO. 
arguments: 

ERO  real  *8 
EVO  real  *8 
MDAMP  integer 
EROUT  real  *8 
EVOUT  real  *8 


EVO.  MDAMP.  EROUT.  EVOUT) 

basic  position  explicitness  parameter  eron  entry 

basic  velocity  explicitness  parameter  e , on  entry 

number  of  cycles  +1  over  which  initial  filtering  by  ADINC  is  desired 

on  exit  the  most  recent  value  of  er  used  by  ADINC 

on  exit  the  most  recent  value  of  e,  used  by  ADINC 


The  next  section  discusses  the  use  of  SETEPS  to  vary  the  explicitness  parameters  The 
arguments  EROUT  and  EVOUT  allow  the  user  to  query  ADINC  about  the  current  values  of 
«,  «,  being  used  as  well  as  to  change  them.  The  need  to  monitor  the  initial  filtering  feature 
(discussed  in  the  next  section)  and  the  need  to  print  out  the  time  averaged  velocities  and 
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pressures  actually  used  to  advance  the  positions  (Eq.  (5'))  and  the  velocities  (Eq.  (6  ))  makes 
the  inclusion  of  EROUT  and  EVOUT  as  arguments  to  SETEPS  a necessity. 

Entry  ADINCO  is  provided  so  the  user  can  monitor  the  performance  of  ADINC.  Each 
time  ADINCO  is  called,  the  number  of  calls  to  ADINC,  the  number  of  timesleps  performed  by 
ADINC,  and  the  number  of  interations  used  in  these  timesleps  is  printed  out.  Each  call  to 
ADINCO  with  MODE  = 0 reinitializes  these  counters  so  we  only  get  the  accumulated  calls, 
timsteps,  and  iterations  since  the  previous  call  to  ADINCO.  MODE  = 1 controls  a special  ini- 
tial filtering  facility  described  in  the  next  section  and  MODE  > 1 is  currently  not  permitted. 

The  major  routine  ADINC  advances  the  physical  variables  {/•,),  (F,),  (p ,),  and  ( /*, } from 
one  timestep  to  the  next  ADINC  uses  the  geometry  and  equation  of  state  routines  and  con- 
trols all  the  logic  of  the  iteration  and  timestep  subcycling.  The  whole  package  has  been  written, 
as  nearly  as  possible,  to  conform  with  styles  and  usage  of  coupled  ordinary  differential  equation 
packages.  Thus  the  user  can  ask  the  package  to  integrate  over  an  interval  requiring  many 
timesteps.  The  user  can  also  change  the  convergence  criteria  and  the  order  of  the  algorithm  via 
the  e,  and  ev  coefficients.  Unlike  ordinary  differential  equations,  however,  partial  differential 
equations  require  boundary  as  well  as  initial  conditions.  Furthermore,  the  generally  large 
number  of  cells  and  concommitant  large  number  of  coupled  equations  force  corresponding  res- 
traints as  the  algorithm.  Because  the  number  of  simultaneous  equations  is  large,  the  algorithm 
must  be  simple  and  optimized  so  the  computation  time  is  acceptable.  ADINC  is  second-order 
accurate  at  best  but  is  fully  vectorized  for  parallel  computation  on  machines  such  as  the  NRL 
Texas  Instruments  ASC. 

The  algorithm  described  in  Section  III  is  a single-step  algorithm  so  numerous  copies  of 
the  physical  variables  at  different  time  levels  do  not  need  to  be  maintained.  Multiple  level 
predictor-corrector  schemes  are  the  norm  for  ODE  packages  but  require  an  inordinate  amount 
of  computer  storage  when  many  equations  are  coupled  as  in  fluid  dynamic  problems.  The  next 
section  considers  the  use  of  ADINC  for  the  rather  general  case  provided  by  the  test  program  in 
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Appendix  B and  supported  by  the  utilities  there.  These  utilities  perform  auxiliary  functions 
useful  in  the  ADINC  computation  but  not  crucial  to  it  such  as  evaluating  internal  and  kinetic- 
energies  and  initializing  a rather  general  multilayer  test  problem.  These  routines  and  their  used 
are  described  via  comments  in  the  listings  in  Appendix  B.  They  would  not  necessarily  be  used 
when  ADINC  is  being  applied  as  a package  within  the  context  of  a more  complex  computation. 


V.  USING  ADINC 

Appendix  B contains  the  listing  of  a rather  general  and  useful  test  program  for  the 
ADINC  package  as  well  as  the  various  utility  programs  not  called  by  the  user  directly.  Use  of 
most  of  the  ADINC  facilities  is  illustrated  there  and  a thorough  understanding  of  this  program 
and  the  examples  provided  in  Appendices  C,  D,  and  E amounts  to  a thorough  understanding  of 
how  to  use  the  ADINC  package.  Here  I present  a general  discussion  designed  to  give  the 
potential  user  an  overall  view  of  the  calculation  and  such  specific  information  as  does  not 
appear  elsewhere  in  this  write-up.  The  three  test  problems  and  variations  discussed  in  the  next 
section  will  illustrate  better  than  any  general  discussion  how  physical  problems  are  solved  accu- 
rately and  reliably  using  ADINC. 

ADINC  is  applied  to  a user  defined  fluid  problem  via  a series  of  structured  calls  to  the  9 
entries  in  the  three  major  routines  dealing  with  geometry,  equation  of  state,  and  fluid  dynamics. 
Figures  3 and  4 outline  the  correct  sequence  for  this  structured  series  of  calls.  Figure  3 con- 
tains the  computations  and  routine  references  required  to  initialize  the  package  and  Figure  4 
contains  those  calculations  and  calls  to  be  performed  during  the  timestep  loop. 

The  first  step  in  the  initialization  (Fig.  (3))  is  to  establish  the  geometry  within  which  the 
calculations  is  to  proceed.  The  integer  a is  set  and  the  values  (G ) for  / = 1,2, 3, 4, 5 are  defined 
if  required.  Then  entry  SETGEO  is  called  to  transmit  these  values  to  ADINC.  Then  the  initial 
cell  interface  positions  are  set  by  the  user  so  that  a call  to  USEGEO  can  calculate  the  interface 
areas,  cell  center  positions,  and  cell  volumes  for  the  initial  grid  geometry. 

The  third  step  is  to  set  the  desired  explicitness  parameters  er  and  e,.  into  the  ADINC 
package.  This  is  done  via  the  entry  SETEPS  where  the  number  of  timesteps  in  the  initial 
damping  transicent  is  also  communicated  to  ADINC. 

Next  the  physical  problem  on  the  grid  must  be  initialized  and  communicated  to  ADINC 
via  the  variables  in  common  block  /ADICOM/.  The  physical  system  has  already  been  divided 
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Figure  3.  To  Initialize  the  AOINC  Package 


• Determine  geometry  (■•  1,  2,  3,  or  4 and  : Gj ; if  needed) 

CALL  SETGEO  (ALPHA,  GEOMCO) 

• Determine  the  initial  interface  positions  : rj ; and  velocities  ; Vj ; 

CALL  USEGEO  (RAD,  AREA,  RADC,  LAMC,  N) 

to  initialize  the  interface  areas,  cell  centers,  and  volumes 

• Set  the  desired  explicitness  parameters  and  damping  period 
CALL  SETEPS  (ERO,  EVO,  MDAMP,  EROUT.  EVOUT) 

to  pass  these  values  on  to  ADINC 

• Determine  the  desired  fluid  properties  and  constants  which  are  conserved 
with  the  cell  motion 

CALL  SETMAT  (MATER,  MASS.  GAMMA,  RHOCON,  N) 
to  set  these  cell  constants  into  common  block 
/ADICOM/ 

• Determine  the  initial  cell  densities  and  pressures 
CALL  SETEOS  (RHO.  PRE,  N) 

to  set  the  current  value  of  the  cell  entropies  ENTC  into  /ADICOM/.  This  has 
to  be  repeated  in  the  timestep  loop  whenever  irreversible  and  dissipative 
phenomena  change  the  entropy. 

F igure  3 — Schematic  lor  initialising  the  ADINC'  package  The  rather  general  case  presented  mirrors  the  test  program 
of  Appendix  B The  test  program  presents,  however,  a number  of  additional  facilities  lor  use  in  a scientific  research 
program  which  are  not  part  of  ADINC  proper  These  do  not  appear  in  the  schematic 
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Figure  4.  Application  of  ADINC  During  a Timestep 


• Establish  and  limit  a monotonically  increasing  cycle  number  — 

DO  9999  ISTEP  = 1,  MAXSTP  (in  the  test  program) 

• Set  a timestep  DELTAT  which  depends  on  geometry,  (low  and  other 
considerations.  This  involves: 

CALL  DTFLOW  (RAD,  DABSV,  DTVAL,  N) 
where  DABSV  is  an  array  containing  the  maximum  absolute  value  of  the 
velocity  in  the  previous  two  timesteps  (DABSV  VEL)  works  fine  in  most 
circumstances.  Other  timestep  limiting  calculations  should  also  be 
performed  at  this  time 

• Perform  diagnostics  on  the  geometry  and  physical  variables.  This  usually 
involves  another 

CALL  USEGEO  (RAD,  AREA,  RADC,  LAMC,  N) 
to  update  AREA,  RADC,  and  LAMC  consistent  with  current  cell  interface 
locations  in  RAD 

any  other  diagnostics;  Including  ERGPRT,  are  performed 

• Advance  the  hydrodynamic  variables  by  first: 

setting  the  new  boundary  positions  RLNEW,  RRNEW  and  velocities 
VLNEW,  VRNEW  before 

CALL  ADINC  (RAD,  VEL,  RHO,  PRE.  N,  DELTAT,  ISTEP.  RRNEW. 

RLNEW,  VRNEW,  VLNEW) 

• Reset  the  6r,  6v  values  by  changing  EPSRO.  EPSVO  in 
CALL  SETEPS  (EPSRO,  EPSVO,  MDAMP,  EPSR.  EPSV) 

9999  CONTINUE  (go  back  for  another  timestep) 

Figure  4 — Schematic  for  using  the  ADINC  package  to  integrate  the  Lagrangian  equations  of  motion  The 
upper  portion  of  the  figure  shows  the  use  of  ADINC  in  diagnostics  and  the  lower  portion,  its  use  in  the 
integration  proper  Additional  physics  can  be  added  via  timestep  splitting  This  is  indicated  in  the  test  pro- 
gram 


into  a number  of  computational  cells.  For  each  of  these  Lagrangian  cells  a number  of  quant 
ties  must  be  established  which  are  constant  throughout  the  ADINC  integration.  These  can  be 
initialized  directly  into  the  /ADICOM/  variables  MATERC,  MASSC,  GAMMaC,  and  RHOC 
by  the  user  or  auxiliary  variables  MATER,  MASS,  GAMMA,  RHOCON  can  be  setup  and  the 
entry  SETMAT  called  to  place  these  quantities  into  the  common  block  in  the  desired  locations 
The  test  program  fills  /ADICOM/  directly  but  demonstrates  the  call  to  SETMAT  anyway. 

Remember  that  the  AM-1  = NCELLS+1  interfaces  and  associated  interface  quantities 
reside  in  the  first  through  the  AM- 1-st  locations  of  interface  arrays  but  the  ,V  cells  and  associ- 
ated cell  quantities  must  reside  in  locations  2 through  AM-1.  This  offset  leaves  "cell"  1 and 
"cell"  AM- 1 available  for  fancy  boundary  conditions  or  other  ghost  cell  applications.  In  the  test 
problem,  printouts  of  Appendices  C,  D and  E these  undefined  quantities  appear  as  the  number 
2.  x lO68  or  I1IIII  or  ******  depending  on  the  print  out  format.  Leaving  these  undefined  quan- 
tities in  the  arrays  has  no  adverse  effects  in  the  calculations,  as  will  be  seen,  and  does  have  a 
potential  advantage  Any  off-by-one  errors  introduced  in  new  code  by  the  user  will  show  up  as 
arithmetic  exceptions  involving  the  use  of  these  undefined  cell  quantities. 

Once  the  initial  density  and  pressure  in  each  cell  have  been  set.  the  cell  entropy  ENTC 

$ 

can  be  determined  and  placed  in  /ADICOM/.  The  entry  SETEOS  performs  this  calculation 
using  the  equation  of  state  function.  A separate  entry  is  provided  here  for  efficiency  because  1 
anticipate  resetting  the  entropy  whenever  non-ideal,  dissipative,  or  energy  source  terms  are 
present.  In  more  complicated  systems  such  as  detailed  chemically  reactive  flows  even  the 
values  of  (y  ) (GAMMAC)  change.  The  utility  DUBLOG  is  provided  and  used  in  SETEOS  and 
USEEOS  because  the  current  ASC  sysem  has  an  error  in  the  double  precision  logarithm  routine 
which  prevents  vectorization.  Efficiency  in  the  SETEOS  and  USEEOS  calculations,  particularly 
the  latter,  is  extremely  important  because  these  transcendental  calculations  seem  to  dominate 
the  ADINC  execution  time. 

The  repetitive  computations  within  an  ADINC  simulation  cycle  are  listed  in  Fig.  4 The 
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first  part  of  a cycle  is  to  find  the  correct  timestep,  8r  (DELTAT  in  the  test  program).  Entry 
DTFLOW  is  provided  in  the  ADINC  package  to  calculate  an  estimated  timestep  which  will 
prevent  interface  crossing  but  several  other  timestep  limitations,  which  are  normally  problem 
dependent,  have  to  be  included.  Even  though  the  sound  speed  is  nominally  not  part  of  the 
timestep  calculation,  it  certainly  affects  the  accuracy  and  the  tendency  toward  non-linear  insta- 
bility. The  test  program  illustrates  a number  of  the  considerations  involved  in  timestep  selec- 
tion. 

Once  the  timestep  has  been  chosen  1 usually  recommend  that  all  I/O,  dumps,  and  diag- 
nostics be  performed.  Since  actual  integration  of  the  first  timestep  has  not  ocurred,  performing 
the  diagnostics  here  within  the  timestep  loop  allows  the  initial  conditions  to  be  printed  when 
ISTEP  = 1.  Thus  diagnostic  tests  and  subsidiary  calculations  have  to  appear  only  one  place  in 
the  code.  As  part  to  these  diagnostics  the  geometry  variables  should  be  updated  as  well  via 
USEGEO 

Before  calling  ADINC  to  integrate  the  fluid  dynamic  equations  in  the  loop,  the  boundary 
conditions  for  the  integration  step  must  be  established.  The  values  of  RRNEW,  RLNEW, 
VRNEW,  and  VLNEW  convey  this  information  to  ADINC.  On  entry  {/•(,  (V),  Ip  ),  and  (P] 
contain  the  physical  quantities  at  time  t , the  beginning  of  the  integration  step  ADINC  is  being 
asked  to  integrate  an  internal  8r  up  to  time  r + 8f.  At  this  time  the  interior  values  of  (r|  and 
( V,)  will  be  determined  but  the  boundary  velocities  T,  and  Tv+1  and  hence  the  boundary  posi- 
tions depend  on  external  factors.  Thus 
| RRNEW  = rn+l  ( t + 8/), 

RLNEW  = r,  (/  + 8r), 

VRNEW  = Vn+]  (/  + 8/), 

VLNEW  = V,(t  + 8f ). 

Nominally  VRNEW  and  VLNEW  can  be  anything  but  for  internal  consistency  RRNEW  should 
f satisfy 
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rv+1(/  + Sr)  = rv+|(/)  + «r8/Kw+l(f)  + (1  — e,)8r  + + 8/) 

with  a similar  equation  for  rt(r  + St).  The  test  program  uses  a simplified  version  of  this  to 
avoid  the  complications  of  calling  SETEPS  again  just  to  set  the  boundary  condition.  Instead  the 
previous  values  «,  and  f,  are  used  as  these  will  apply  unchanged  throughout  the  calculation 
except  possibly  in  the  first  few  cycles. 

If  the  values  (P,)  and  I V)  are  stored  in  auxiliary  arrays  ( V“]  and  { /* "}  before  ADINC  is 
called,  the  average  velocities  and  pressures  used  to  advance  the  interface  positions  and  veloci- 
ties can  be  calculated  as  shown  in  the  test  program.  Even  though  this  information  can  be 
reconstructed  without  special  calculation  if  two  sucessive  timesteps  are  printed  out,  these  inter- 
mediate values  are  the  basis  of  the  Lagrangian  dynamics  so  the  prospective  user  of  ADINC 
should  be  familiar  with  the  contents  and  purpose  of  DO  LOOP  110. 

The  user  of  ADINC  also  has  control  over  the  values  of  the  explicitness  parameters  t,  and 
«,  for  the  interface  position  and  velocity  integrations  When  e,  or  = 0,  the  corresponding 
equations  is  fully  forward  differenced  and  therefore  implicit.  When  er  and  t , are  both  equal  to 

y,  the  calculation  is  fully  centered  and  hence  will  be  second-order  accurate  in  time.  The  use 

of  er  > y or  «„  > y can  lead  to  numerical  instability  and  is  therefore  discouraged  The  effect 

of  using  different  values  for  e, , e,  is  discussed  in  the  next  section  as  part  of  the  sound  wave 
test  problem.  At  any  point  during  the  calculation  the  values  of  er,  e,  can  be  changed  by  calling 
the  entry  SETEPS  with  the  desired  values  as  the  first  two  arguments.  The  default  values  are 
«,  = 0.45,  e,  = 0.45  so  the  calculation  will  be  slightly  damped.  This  setting  is  good  for 
moderately  long  timesteps  relative  to  sonic  transit  times  but  will  eventually  damp  finite  speed 

sound  waves  appreciably  so  the  flexibility  to  change  c,  closer  to  y is  important. 

A special  facility  has  been  built  into  ADINC  to  damp  out  initial  transients  when  long 

timesteps  are  being  used.  The  facility  is  controlled  by  the  third  argument  of  SETEPS,  the 

integer  MDAMP.  When  MDAMP  = 1,  the  facility  is  disabled  making  the  values  of  t, 
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constant  until  the  next  time  they  are  changed  under  user  control  by  calling  SETEPS  with  new 
values. 


Before  describing  exactly  what  this  facility  does,  it  is  important  to  understand  why  such  a 
facility  can  be  necessary  in  some  problems.  Consider  a relatively  slow  flow  in  which  we  wish  'o 

Keep  5/  large  for  computational  efficiency  but  where  t,,  t,  — * is  desired  for  numerical  accu 

racy,  i.e.  we  want  to  keep  the  numerical  damping  small  In  such  a problem  any  initial  condition 
of  pressures,  positions,  and  velocities  can  be  resolved  in  high  frequency  sound  waves  and  slow 
flow  components  such  as  might  be  driven  by  thermal  conduction,  chemical  energy  release,  et ' 
With  e — 1/2,  the  high  frequency  components  will  oscillate  rapidly  and  can  mask  the  desired 
slowly  varying  solution  even  though  these  oscillations  are  stable.  The  problem  arises  because 
the  various  initial  conditions  are  slightly  incompatible  with  each  other  and  this  incompatibility 
does  not  decay  away  quickly  when  e. , e are  near  1/2.  Rather  than  expending  a great  deal  of 
effort  to  filter  the  initial  conditions,  it  is  often  perfectly  adequate  to  damp  the  first  few  cycles  of 
the  ADINC  calculation  strongly  by  using  e,.  e < 


The  algorithm  implemented  linearly  increases  er.  e,  from  ERO/MDAMP,  EVO/MDAMP 
to  ERO.  EVO  in  MDAMP-i  successive  calls  to  ADINC.  The  default  values  are  MDAMP  --- 
10.  ERO  = 0.45,  EVO  = 0.45  if  the  SETEPS  entry  is  not  used.  The  implementation  i>  ,.ie 
obvious  in  Appendix  A.  Clearly  MDAMP  - 1 circumvents  the  changing  values  ol  t.,  t. 
entirely.  Since  it  may  be  desirable  to  perform  this  filtering  operation  more  than  once  during  a 
calculation,  the  facility  is  provided  to  reset  the  counter  for  the  MDAMP-1  filtering  steps.  Every 
t-  v \DINCO  is  called  with  .MODE  = 0,  the  integer  NCALL  gets  reset  to  zero  after  the  three 
integers  are  printed  out.  When  MODE  > 1,  NCALL  does  not  get  reset  and  hence  the  filtering 
operation  is  not  performed  anew. 


The  ADINCO  entry  was  provided  so  the  user  can  diagnose  the  performance  of  ADINC  at 
any  point  in  a calculation.  When  called,  the  total  number  of  calls  to  ADINC,  the  total  number 
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of  subcycle  steps  performed,  and  the  total  number  of  iterations  performed  since  the  last  called 
to  AD1NCO  are  printed  out  Each  time  ADINCO  is  called  with  MODE  — 0,  the  counters  are 
all  reset.  Using  ADINCO  to  control  re-initializaiion  of  the  filtering  facility  MDAMP  ensures 
that  diagnostic  prints  are  obtained  every  time  the  solution  is  re-filtered. 
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VI.  DISCUSSION  OF  TEST  PROBLEMS 


Appendices  C,  D,  and  E provide  selected  printouts  from  three  test  calculations  using 
ADINC.  These  test  are  chosen  to  illustrate  the  use  of  ADINC  in  several  different  and  interest- 
ing cases,  to  demonstrate  the  flexibility  and  accuracy  of  the  package,  to  demonstrate  some  of 
the  control  features  the  user  has  at  his  command,  and  to  allow  positive  code  * erification  via 
known  answers  for  users  implementing  the  code  for  themselves.  In  addition  to  providing  evi- 
dence that  ADINC  and  its  associated  subroutines  do  something  reasonable  and  accurate  in 
several  cases,  these  test  problems  provide  a jumping  off  point  for  many  user  applications  and 
may  be  directly  adaptable  with  little  or  no  reprogramming. 

The  three  problems  are: 

#1.  An  Adiabatic  Sound  Wave  Test, 

#2  An  Incompressible  Slug  Between  Adiabatic  Gases,  and 

#3.  A LINUS  Simulation. 

Each  tes'  has  a main  calculation,  for  which  results  are  printed  in  the  appendices,  along  with  the 
corresponding  data  for  code  verification.  In  this  section  a number  of  auxiliary  calculations  are 
also  reported  to  establish  properties  of  the  ADINC  package  with  respect  to  convergence,  accu- 
racy, and  stability. 

#1.  An  Adiabatic  Sound  Wave  Test  (see  Appendix  C) 

The  first  test  problem  concerns  the  ability  of  ADINC  to  capture  the  properties  of  a sound 
wave  in  an  ideal  gas.  The  gas  is  uniform  with  density  1.4  and  pressure  1.0  in  the  absence  of 
the  sound  wave.  The  gas  constant  y = 1.4  ensures  a sound  speed  C,  = 1.0=  \TyPjp  ■ The 
gas  is  contained  between  two  rigid,  impermeable  walls  at  ,v=0  and  v = 1.0  so  a sound  wave  pro- 
pagates from  x=0  to  x = l and  back  in  a period  of  time  r=2.0.  The  standard  test  #1  has  hall  a 
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wavelength  between  the  walls,  the  velocity  being  initialized  sinusoidally  according  to 

V(x.o)  — 8 V sin  nx  (28) 

where  bV(DVEL)  = 0.01  for  the  standard  test.  The  initial  grid  is  uniform  with  8.v  = 0.1  and 

the  timestep  is  8 z =0.1.  The  standard  calculation  is  thus  performed  with  20  cells  per 

wavelength  and  20  timesteps  per  period. 


The  data  for  this  standard  test  #1  are  given  in  the  following  table: 


Table  1.  ADINC  Input  Data  for  Standard  Test  #1  — Adiabatic  Sound  Wave 


Namelist  /CONTRL/ 
MAXSTP  = 26 
IPRINT  = I 
ALPHA  = 1 
N = 10 


(program  control) 

DTMIN  = 0.1 
DTMAX  = 0.1 
EPSRO  = 0.5 
EPSV0  = 0.5 
GEOMCO  = (1.0,0.0,0.0,0.0,0  0)  (not  used  if  <r  = l) 
LZONE  = FALSE.  LTCND  = FALSE. 
LCHEM  = FALSE.  LTPRT  = TRUE. 

LD1FF  = FALSE. 


Namelist  /SHL1NI/ 
NSHELL  = 1 
RN  = 1.0x10 
VN  = 0.0 


(layer  (shell)  initializer) 
MODE  = I 
DRHO  = 0.0 
DVEL  = 0.01 


Namelist  /SHLDAT/ 
LCELLS  = 10 
RN  = 1.0 
VN  = 0.0 
RHOS  = 1.4 
POWS  - 1.0 


(first  and  only  layer) 
MATERS  - 1.01 
GAMMAS  = 1.4 
RHOCS  = 0.0 
PRES  = 1.0 


Because  the  peak  velocity  is  small  compared  to  the  sound  speed,  8 V/C\  =10  \ the  evolu 
tion  of  nonlinear  effects  will  be  slow  and  the  profiles  will  remain  sinusoidal  for  a long  tune 
The  period  of  the  oscillation  is  determined  by  noting  when  rh,  which  starts  at  0.5,  passes 
through  0.5  the  second  time.  Since  the  numerical  period  computed  by  ADINC  will  differ  from 
the  exact  theoretical  period  due  to  finite  difference  truncation  error,  interpolation  is  necessary 
between  timesteps  to  find  the  numerical  period.  This  interpolation  process  has  an  error  associ- 
ated with  it  so  the  numerical  period  is  never  found  exactly. 

The  standard  test  #1  with  20  cells  per  wavelength  and —20  timesteps  per  cycle  has  the 
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crossmi;  ol  rh  through  0.5  occur  between  cycles  20  and  21  At  r-2.0.  rk  - 0 4991b  and  at 


| 

L 

F 


f— 2.1  the  code  gives  />  — 0.50074.  The  period  computed  Irom  these  values  is  ra>dc  —■ 
2.0245  ± .001,  in  error  by  ~ 1%.  This  error  arises  Irom  several  sources,  the  most  important  of 
which  are  the  finite  timestep  and  grid  size  of  the  calculation.  The  error  from  interpolating  for 
r„=0  5 is  twenty  times  smaller. 

1 able  2 below  summarizes  the  period  computed  using  ADINC  for  several  values  of  grid 
size  (initial)  and  timestep.  The  number  of  cells  per  wavelength  and  limesteps  per  period,  the 
important  non-dimensional  quantities,  are  also  shown 


Table  2.  Numerically  Computed  Period  for  Adiabatic  Sound  Wave 


Since  the  values  of  t,  and  e,.  are  both  0.5,  the  time  integration  is  centered  as  well  as  the  spatial 
differences  Therefore  we  expect  full  second  order  accuracy.  The  coefficients  in  time  and  space 
are  differeni  ol  course  since  the  spatial  and  temporal  algorithms  are  different.  The  functional 
form  of  the  computed  period  is 

^ ^ ^ ( 1 Tot f fix ) ^ ) ( 1 f S r ) -h  higher  order  terms  (29) 

where  a~-l  6 and  li  — 3 2 when  8x  and  8/  are  measured  in  units  of  wavelengths  and  wave 
periods  respectively.  That  the  above  data  in  Table  2 satisfy  Eq.  (29)  can  be  seen  by  the  fact 
that  the  error  in  the  20  20  calculation  is  four  times  the  error  in  the  40-40  calculation. 

In  these  calculations  the  layer  summaries  printed  by  the  test  program  include  energy 

sums.  Subroutine  ERGPRT  and  the  diagnostics  section  of  the  main  program  show  how  these 

sums  are  calculated.  The  interval  energy  density  Ethr„  is  given  by 
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Eri*m,  -P./lV,' "I)  (30' 

for  each  cell  of  volume  A The  kinetic  energy  density  is  a little  more  difficult  to  compute 
because  the  interface  velocities  are  known,  not  the  cell  center  velocities.  Averaging  the  veloci- 
ties to  the  cell  center  introduces  a damping  term  which  appears  as  an  overall  small  amplitude 

* 

oscillation  in  the  total  energy  as  the  wave  cycles  between  potential  and  kinetic  energy. 

There  is  a definition  of  the  kinetic  energy  which  is  more  consistent  with  ADINC’s  finite 
difference  algorithms  and  which  does  not  display  the  energy  changes  that  the  velocity  average 
definition  does.  This  simple  diagnostic  uses  the  fact  that  ADINC  defines  cell  center  positions 
{/?,)  and  matches  accelerations  across  cell  interfaces.  Let  A ' ana  A be  the  cell  partial 
volumes  to  the  right  and  the  left  of  the  cell  center  position  respectively 

a = v * + a,-  no 

In  the  test  program  |A  +l  and  (A,  ) are  calculated  as  a simple  average  but  more  accurate  values 
can  be  determined  by  USEGEO  for  non-Cartesian  geometries.  The  kinetic  energy  density  in 
each  cell  is  given  by  the  formula 

Ek,ne,  A,  = yp  [k2A+  + Ki,  A ] (32) 

The  code  diagnostics  print  out  the  total  energy  to  nine  significant  digits  and  to  this  accu- 
racy energy  is  conserved  identically  Since  the  algorithm  :s  nominally  reversible,  good  energy 
conservation  is  expected  for  this  sound  wave  test  problem.  Fssent  illy  perfect  energy  conserva- 
tion here  means  consistent  thermal  and  kinetic  energy  definitions  have  been  found. 

A separate  calculation  with  20  cells  was  performed  In  which  the  individual  zones  varied  in 
thickness  by  a factor  of  100,  alternating  between  fix  =•  - and  fix  = yy  all  the  way  irom 

x=0  to  x = l The  timestep  was  taken  to  be  8/  = 0.05  and  the  period  was  determined  to  be 
2.0121,  In  error  by  ~ 0.5%.  The  conservation  of  energy  for  this  calculation  was  again  good  to 
I part  in  10"  even  though  ADINC  had  to  subcycle  twice  at  each  timestep  because  some  of  the 
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cells  were  so  small.  The  wide  disparity  in  cell  sizes  did  not  adversely  affect  the  accuracy  of  the  calcu- 
lation. Of  course  the  large  cells,  not  the  small  ones,  determine  the  overall  accuracy  of  the  cal- 
culation but  more  realistic  problems  with  widely,  varying  characteristic  scale  lengths  and/or 
imbedded  discontinuities  will  be  able  to  take  good  advantage  of  this  additional  flexibility  in  the 
ADINC  algorithm. 

A series  of  calculations  was  next  performed  (with  uniform  zoning  again)  where  the  values 
of  (,  and  t,  were  varied  to  determine  the  amount  of  sound  wave  damping  implied  by  forward 
differencing.  It  is  important  to  understand  that  the  stability  properties  usually  invoked  as 
benefits  of  implicit  differencing  come  at  stiff  price.10"13  A user  of  ADINC  should  perform  cal- 
culations with  t,  t , as  close  to  0.5  as  possible  so  accuracy  can  be  maintained  as  well  as  stability. 
To  measure  wave  damping  the  energy  diagnostic  was  used  to  filter  out  the  oscillations.  Twenty 
cells  per  wavelength  and  20  timesteps  per  cycle  were  again  chosen  as  the  standard  conditions. 
Figure  5 and  Table  3 summarize  a number  of  computations  performed  with  different  values  of 
t,  and  the  ADINC  explicitness  parameters.  The  wave  amplitude  relative  to  the  initial  wave 
amplitude,  „ is  plotted  versus  the  time  measured  in  wave  periods  (theoretical).  Jilt)  is 
determined  by  subtracting  from  the  total  energy,  the  thermal  energy  that  the  system  would 
have  with  zero  velocity.  The  linearity  of  versus  t/r  on  the  semi-logarithmic  scale  demon- 
strates that  the  numerical  damping  is  exponential  as  expected.  Even  with  the  very  modest 
damping  introduced  by  — 0.45,  the  sound  wave  has  decayed  by  more  than  a factor  of 

2.5  in  only  10  oscillation  periods.  This  means  that  the  use  of  fully  forward  differenced  schemes 
for  numerical  stability  makes  the  meaningful  simulation  of  sound  waves  essentially  impossible. 

Table  3 below  summarizes  these  damping  calculations  in  terms  of  the  per  period  damping 
coefficient.  The  physical  sound  wave  being  studied  should  undergo  no  damping  at  all  so  the  loss 
of  energy  marks  a numerical  error  which  is  used  here  to  measure  the  accuracy  of  the  ADINC 
algorithm  Variation  of  damping  with  t,  timestep  8 1,  and  grid  size  8c  is  shown. 
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(t/r)  -* 

Figure  5 — Wave  amplitude  as  a function  of  time  illustrates  damping  introduced  by  forward  differencing  Amplitudes 
are  measured  in  units  of  the  initial  amplitude  The  values  of  the  explicitness  parameters  are  shown.  With  t,  = r,  = 
0 45.  for  example,  the  wave  has  damped  by  a factor  of  >2.5  in  ten  periods  of  the  wave. 
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Table  3.  Variation  of  wave  damping  coefficient  (per  period)  wilh  numerical 
control  parameters.  The  physical  wave  should  be  undamped 


cells  per 
wavelength 

timesteps 
per  period 

damping 

coefficient 



ESM 

0.50 

mm i 

■KiHH 

~ HI' 

mm 

0.475 

JrH 

wm 

mm 

0.45 

20 

20 

.909 

wrm 

0.40 

20 

20 

.826 

0 30 

0.30 

20 

20 

.683 

0.40 

0 SO 

20 

20 

K3 

0.50 

0.40 

20 

20 

mm 

0.40 

0.40 

10 

20 

.830 

0.40 

0 40 

20 

40 

— ■ 

There  are  no  surprizes  in  these  data.  They  are  a quantitative  as  well  as  qualitative  verification 
of  the  linearized  ADINC  dispersion  relation 
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(34) 


These  formulae  predict  the  symmetric  behavior  of  wave  damping  and  wave  period  with  e,  and 
<,  they  also  predict  the  relative  insensitivity  of  the  wave  damping  to  the  cell  size  when  8a-  is 
small  in 


_ 2(7,8/  *8x 

/3  = — — sin  — - 
8a  2 


(35) 


The  next  set  of  tests  performed  concerned  the  nonlinear  behavior  of  sound  waves  com- 
puted by  ADINC.  Three  variations  of  the  standard  test  #1  were  performed  in  which  8)',  the 
initial  sinusoidal  velocity  maximum,  was  increased  to  0.1,  0.3.  and  0 5 limes  the  speed  of  sound 
in  the  system.  Figure  6 shows  four  normalized  velocity  profiles  at  cycle  40  in  the  three  calcula- 
tions, two  linear  wave  periods  after  the  calculation  is  initialized.  The  linear  (sinusoidal)  profile 
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Figure  6 — Nonlinearity  in  Standard  Test  #1  Without  damping,  increasingly 
nonlinear  initial  sound  waves  develop  short  wavelength  oscillations  propor- 
tionately faster  The  profiles  are  scaled,  the  sinusoidal  linear  solution  after  two 
periods  is  also  shown  for  comparison 


is  shown  as  the  solid  line  for  comparison. 


In  ihe  ft!  =0.1  case  the  calculation  required  more  iterations  than  in  the  previous  cases 
and  ADINC  had  difficulty  with  convergence  criteria  during  some  of  the  timesteps.  The  choice 
of  t,  = t,  = 0 5 for  these  calculations  seems  to  be  permitting  the  nonlinear  growth  of  short 
wavelength  perturbations  which  are  not  preferentially  damped  in  the  calculation.  At  propor- 
tionately larger  wave  amplitudes  these  perturbations  get  larger.  In  the  &F  = 0.5  undamped 
case,  strong  evidence  of  an  even-odd  type  of  grid  instability  can  be  seen.  Clearly  some  damp- 
ing is  necessary  for  nonlinear  cases  like  these  even  before  shocks  can  be  expected  to  form. 

The  Sf  = 0.5  calculation  was  repeated  with  0.45,  with  30  cells  in  the  system 

(60  cells  per  wavelength),  and  with  Sr  = 0.02  Figure  7 shows  the  velocity  profile  at  a series  of 
times  during  the  calculation  half  a theoretical  linear  wave  period  apart.  In  the  physical  problem 
simulated  here  two  large  amplitude  sound  waves  of  equal  amplitude  are  passing  through  each 
other  travelling  in  opposite  directions.  As  the  pressure  increases  in  the  wave  so  does  the  sound 
speed.  Thus  the  two  waves  each  steepen  until  a shock  (actually  a number  of  shocks)  forms. 
Even  though  there  is  some  linear  damping  in  the  calculation  of  Fig.  7,  it  is  not  enough  to  pro- 
vide the  required  conversion  of  kinetic  to  thermal  energy  through  the  shock.  Future  emphasis 
of  research  on  the  ADINC  package  itself  will  be  on  these  nonlinear  aspects  of  implicit  hydro- 
dynamics. This  is  a very  difficult  problem.  Implicit  differencing  means  that  the  numerical 
"characteristic"  travels  at  infinite  speed.  Thus  the  possibility  always  exists  for  information  about 
the  shocked  fluid  to  propagate  nonphysicallv  ahead  of  the  advancing  shock.  My  current  recom- 
mendation for  treating  shocks  accurately  is  to  use  the  FCT  algorithms  described  in  refs.  14-16. 

#2.  An  Incompressible  Slug  Between  Adribatic  Gases  (see  Appendix  D) 

The  second  test  problem  is  illustrated  in  Figure  8.  There  are  three  layers  bounded  on  the 
left  and  on  the  right  by  rigid  impermeable  walls.  The  center  region  is  an  incompressible  slug. 
The  left  and  right  regions  are  each  adiabatic  gas  layers  of  exactly  the  same  properties  as  the  sin- 
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NONLINEAR  SOUND 
WAVE  STEEPENING 


- t = 0.0  (LINEAR  WAVE) 

• t = 1.0  = t/2 
■ t = 2.0  = t 

* t = 3.0  = 3t/2 


Figure  7 — Velocity  profiles  of  nonlinear  sound  wave  steepening  The  four  curves  are  taken  half  a 
(theoretical  linear)  period  apart.  The  r/2  and  3t/2  data  are  inverted  + for  - and  flipped  left  for  right  on 
the  plot  so  that  the  nonlinear  evolution  is  evident  With  t,  — e,  - 0 45  as  in  this  calculation,  there  is 
not  enough  dissipation  to  maintain  monotonicity  at  the  shock 


gle  gas  layer  considered  as  standard  test  #1.  The  large  density  discontinuities  across  the  two 
layer  interfaces  makes  this  a test  of  the  acceleration-matching  part  of  the  ADINC  algorithm. 
Without  this  correction  the  effective  mass  of  the  slug  would  vary  by  half  a cell  or  so  hence 
finite  difference  errors  would  be  of  order  5%. 

Detailed  comparisons  require  an  analytic  solution.  Let  the  slug  have  a mass  per  cm' 
and  a density  ps.  The  right  and  left  hand  adiabatic  gas  layers  have  densities  which  are  deter- 
mined by  the  location  of  the  moveable  slug.  At  the  initial  time  t = 0 the  slug  is  at  the  equili- 
brium position  but  moving  with  a nonzero  velocity  Vs.  The  initial  pressures  are  equal, 
PR  = P(,  but  the  densities,  pR  and  pL , the  width  of  the  regions,  WR  and  WL,  and  the  gas  con- 
stants, yR  and  yL , can  all  be  different.  This  system  is  soluble  both  in  the  linear  and  in  the 
asymptotic  nonlinear  limit  and  thus  permits  a detailed  evaluation  of  the  accuracy  of  the  numeri- 
cal techniques. 


Let  the  displacement  of  the  slug  from  equilibrium  be  denoted  by  X(t).  Then 


WLU)  « WL+XU)  PLU)  = SLWL^ L(z) 


(36) 


WRU)  = WR-XU),  PRU)  = SRWRr*(t) 

where  the  fixed  entropies  SL  and  SR  can  be  determined  from  the  initial  conditions.  There  are 


two  equations  to  be  solved. 


dXU ) 
dt 

dVU ) 
dt 


V(t)  and 

PLU)  - PRU). 


(37) 


Linearizing  gives  PL(t ) ~ Pl  ~ y lsl^'l  l ' X.  With  a similar  expression  for  PRU).  The 
equilibrium  terms  cancel  and 


M , » -XU) 

dt 


yt.Pi.  yRpK 


Wn 


(38) 


(39) 


1 
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The  test  calculations  performed  on  this  problem  with  ADINC  have  W w = W,  = 1.0cm 
and  the  slug  density  p ,=  140  gm/cm',  the  slug  length  is  1.0  cm,  and  M.  - 140  gm/cm  The 
initial  gas  layer  pressures  are  1.0  dyne/errr  and  yL  = y^  = 1-4.  With  these  parameters  the  fre- 


quency co  = v'2. 8/140  = 141421  and  the  period  t = = 44.4288  is  expected. 

OJ 


Actually  we  should  not  neglect  the  mass  of  the  gas  because  ADINC  won’t.  When  the 
density  ratio  is  large  as  it  is  here,  the  velocity  will  vary  roughly  linearly  across  the  gas  regions  in 
order  that  the  density  and  pressure  stay  spatially  constant.  The  elective  mass  of  the  slug  there- 
fore increases  by  half  of  the  mass  of  each  ot  the  gas  regions.  Equation  (39)  becomes 

^ = lyL  Pl/Wl  Pr/Wr)  (40) 

[Ms  + — M/  + — Mr] 

which  predicts  a period  r = 44.6504  sec.  Here  I will  not  attempt  still  higher  order  corrections 
which  arise  from  the  curvature  of  the  gas  layer  profiles,  i.e.  the  fact  that  the  density  cannot  be 
constant  if  there  is  a pressure  variation  across  the  gas  layer.  These  are  of  order  1:10  4 and  are 
about  as  large  as  the  interpolation  error  in  evaluating  the  data,  ±0.001  sec. 


The  data  for  the  standard  test  #2  are  given  in  the  following  table: 


Table  4.  ADINC  Input  Data  for  Standard  Test  #2-Incompressible 
Slug  Between  Adiabatic  Gases 


Namelist/CONTRL/  (program  control) 

MAXSTP  = 26  DIMIN  = 1 .0 

IPR1NT  = 1 DTMAX  = 1.0 

ALPHA  - 1 EPSR0  = 0.5 

N = 10  EPSV0  = 0.5 

GEOMCO  = (1.0,  0.0,  0.0,  0.0,  0.0)  (not  used  if  «=1) 

LZONE  = FALSE.  LTCND  = FALSE. 

LCHEM  = FALSE.  LTPRT  = TRUE. 

LDIFF  = FALSE. 

Namelist  /SHLINI/  (layer  (shell)  initializer) 

NSHELL  - 3 MODE  = 1 

RN  = 1.0xl0~2°  DRHO  = 0.0 

VN  = 0.0  DVEL  - 0.0 

Namelist  /SHLDAT/  (first  layer  data) 

LCELLS  = 5 RN  = 1.0 

MATERS  =1.01  POWS  = 1.0 


VN  = 0.001 
RHOCS  = 0 0 


1 


r- 


GAMMAS  = 1.4 

Namelist  /SHLDAT'  (second  layer  data) 
LCELLS  = 5 
MATERS  = 2.02 
GAMMAS  = 0.5 


Namelist  /SHLDAT/  (third  layer  data) 
LCELLS  = 5 
MATERS  = 3.01 
GAMMAS  = 1.4 


PRES  =1.0  RHOS  — 1.4 


RN  = 2.0  VN  = 0.001 

POWS=1.0  RHOCS  = 140 

PRES  = 10 

RHOS  = 140.000000000140 


RN  = 3.0  VN  = 0.0 

POWS  = 1.0  RHOCS  = 0.0 

PRES  =1.0  RHOS  =1.4 


In  the  second  layer,  absolute  incompressibility  is  not  forced  so  ii  can  be  used  as  a diagnos- 
tic of  ADINC’s  accuracy.  The  density  in  this  slug  satisfies 

p,  = 140  + (P/S)2  (41) 

where  the  initial  value  of  ps  = 1 40 [ 1 + 1 0 1 2]  is  used  in  conjunction  with  pL  = 140  in  the  ini- 
tializer to  determine  S = 8.4516xl04.  In  this  second  test,  therefore,  density  fluctuations  of 
order  1(T16  are  expected  physically.  This  is  smaller  than  roundoff  error  The  calculations 
described  below  all  have  the  slug  density  constant  to  at  least  1 part  in  10ll!,  about  the  conver- 
gence criterion  for  the  nonlinear  iteration  in  ADINC. 


The  ADINC  package,  of  course,  accepts  a wide  continuous  range  of  values  for  p > 
and  S in  the  formulae  for  the  equation  of  state  but  for  only  a few  values  of  gamma  is  the 
energy  integral  readily  computable  (I  believe).  When  p,  = 0 in  Eq.  (3).  tor  example,  the 
energy  density  £is  given  by  the  familiar  formula 

E = P/(y—\ ).  (42' 

when  pi  ^ 0,  £ is  given  by  the  integral 

i 

f=y  J (p(L')  — ot)y  dL  (43) 

where 

p( L ) = p, LJL’  (44) 

and  L,  is  the  size  of  the  system  when  p=p(.  The  integral  (43)  is  basically  just  an  evaluation  of 
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the  f PdV  work  done  on  the  fluid. 


When  y = 1/2,  as  in  the  test  case,  the  integral  can  be  performed  in  closed  form. 

[ja_  tan~'(p/p,  — 1) 1/2 


P, 


(p/p,-D 


1/2 


- 1 


(45) 


This  formula  for  the  internal  energy  density  goes  to  zero  in  the  zero  pressure  p = p,  limit.  As 
will  be  seen  in  the  test  calculation  of  Appendix  D,  the  interna!  energy  calculated  for  the  slab  is 
thirteen  orders  of  magnitude  smaller  than  the  thermal  energy  of  the  gas  layers.  This  indicates 
near  perfect  effective  incompressibility  of  some  finite  difference  cells  immediately  adjacent  to 
cells  whose  equation  of  state  requires  compression.  Equation  (45)  is  implemented  as  an  energy 
diagnostic  for  all  cells  having  y—  1/2  in  the  test  program  so  that  future  users  may  see  examples 
of  how  the  solution  quality  vcries  as  a result  of  varying  both  physical  and  numerical  parameters 
in  the  calculation. 


Again  temporal  and  spatial  resolution  are  expected  to  play  major  roles  in  determining  the 
accuracy  of  the  computed  solution.  Table  5 below  summarizes  the  results  of  a number  of  test 
calculations  with  varying  timestep  and  cell  size. 


Table  5.  Numerically  Computed  Period  for 
Oscillatory  Incompressible  Slug 


As  can  be  seen,  the  timestep  plays  a much  larger  role  here  in  the  accuracy  than  the  spatial  reso- 
lution. The  errors  in  the  temporal  integration  are  about  the  same  as  for  the  pure  sound  wave 
when  5/  is  measured  in  units  of  the  appropriate  period.  Here  the  slug  is  very  heavy  so  the 
period  is  much  longer  than  in  test  #1  reported  above.  At  —20  timesteps  per  period  the  error  is 
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again  between  0.5%  and  1%  even  though  the  period  is  theoretically  more  than  22  limes  longer 
for  the  slug  than  for  the  sound  wave.  This  suggests  a usefull  rule  of  thumb  Expect  < 1 % 
errors  in  computation  phenomena  with  a characteristic  period  of  twenty  finite  difference  umesteps. 
Faster  phenomena  will  be  less  accurate,  slower  phenomena  more  accurate.  Computed  periods 
seem  to  be  slower  titan  theoretically  expected. 

The  acceleration  matching  algorithm  in  ADINC  is  the  reason  that  the  computed  slug 
period  is  insensitive  to  the  spatial  resolution.  The  interfaces  are  perfectly  resolved  and  any  den- 
sity discontinuities  are  automatically  treated  with  high  accuracy  as  long  as  pressure  variations 
are  linear  between  cell  centers  and  adjacent  interfaces. 


As  with  standard  test  #1,  the  nonlinear  limit  of  this  oscillatory  slug  problem  is  interesting. 
When  the  pressure  in  the  gas  layers  is  low  relative  to  the  kinetic  energy  of  the  slug,  the  slug 
travels  from  one  wall  to  the  other  at  essentially  constant  velocity  where  it  then  rebounds  specu- 
larly. Therefore,  to  lowest  order  anyway,  the  period  of  oscillation  will  be 


2(  WL  -l-  WR) 


(46) 


Equation  (46)  gives  only  the  limiting  value,  of  course,  valid  when  the  gas  layer  compresses  to 
vanishing  thickness  during  rebound.  Figure  9 plots  the  slug  frequency  1/r  determined  numeri- 
cally as  a function  of  the  maximum  slug  velocity.  Clearly  the  correct  asymptote  is  being 
approached.  As  the  slug  velocity  approaches  the  sound  speed  of  the  gas  in  the  adiabatic  layers, 
however,  the  flow  becomes  complicated  and  shock  dissipation  (which  ADINC  cannot  handle 
properly)  will  become  important. 


As  in  standard  test  #1,  further  investigations  of  steepening  sound  waves  will  be  deferred 
until  a future  paper.  New  users  of  ADINC  might  find  it  very  instructive  to  boost  the  sound 
speed  (or  lower  the  pressure)  in  the  gas  layers  to  allow  a faster  slug  velocity  without  generating 
a shock  in  the  buffer  layers.  In  principle  the  density  of  these  layers  can  always  be  dropped  far 
enough  with  />  = 1 so  that  any  given  slug  velocity  T,  will  be  subsonic  throughout  the 
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1/T 


OSCILLATION  FREQUENCY 

vs 


PEAK  SLUG  VELOCITY 


'/slug  (0) 

Figure  9 — Oscillation  frequency  versus  peak  slug  velocity  For  Fs  near  zero  the  fre- 
quency is  essentially  constant  at  the  theoretical  value  given  by  Eq.  (40)  For  large  veloci- 
ties the  slug  always  travels  close  to  hut  below  the  peak  velocity  so  the  actual  frequency- 
lies  somewhat  below  the  asymptotic  limiting  curve 


compression  and  rarefaction  phases  of  the  slug  oscillation.  When  shocks  do  develop,  1 recom- 
mend the  FCT  algorithms  discussed  in  refs.  14-16. 


#3.  A LINUS  Simulation  (see  Appendix  E) 


The  problem  which  necessitated  writing  ADINC  in  its  current  form,  more  than  any  other, 
is  the  simulation  of  the  slightly  compressible  liquid  liners  used  in  Linus6"*.  Original  treatment 
of  the  plasma/magnetic  field/liner  system  assumed  the  cylindrical  plasma/magnetic  field  core 
was  being  compressed  by  a constant  density,  incompressible,  ideal  liquid  flowing  radially  inward 
The  compressibility  of  the  liner  actually  can  have  a substantial  effect  on  the  dynamics,  energet- 
ics and  dwell  time  of  the  implosion.  Hence  the  desire  to  simulate  acoustic  phenomena  accu- 
rately for  time  periods  long  compared  to  a sonic  transit  time  became  a necessity.  Fully  forward 
differenced  schemes  and  even  the  quasi-static  iteration  used  originally  for  Linus10  '3  would  not 
perform  this  calculation  well.  Therefore  ADINC  was  developed  with  the  capability  for  central 
as  well  as  forward  differencing  under  programmer  control 


The  geometry  of  the  calculation  is  shown  in  Figure  10  with  the  reference  values  of  the 
coordinates  diplayed.  Five  cylindrical  shells  are  used  with  a 100  atmosphere  pressure  in  the 
driver  gas  plenum.  The  input  data  for  the  standard  test  #3  are  given  in  the  following  table- 


Table  6 ADINC  Input  Data  for  Standard  Test  #3  — A LINUS  Simulation 


Namelist/CONTRL/  (program  control) 
MAXSTP  - 1001 
IPRINT  = 50 
ALPHA  = 2 
N = 30 


DTMIN  =10  ’ 
DTMAX  =10  4 
EPSR0  = 0.45 
EPSV0  = 0 45 


GEOMCO  = (1.0,  0.0,  0.0,  0.0,  0.0)  (not  used  since  a=2) 

LZONE  = FALSE  LTCND  = .FALSE. 

LCHEM  = FALSE.  LTPRT  = TRUE 

LD1FF  = FALSE. 

Namelist  /SHL1NI/  (layer  (shell)  initializer) 

NSHELL  = 5 MODE  = 1 

RN=10m  drho-o.o 

VN  = 0.0  DVEL  = 0.0 


Namelist  /SHLDAT/  (plasma  layer  data) 
LCELLS  = 3 


RN  = 10.0 


VN  = 0 0 


■ 


1 


r- 


i 


f 

f 


MAI  ERS  - 1 .01 

POWS  - 1 0 

GAMMAS  = 5/3 
PRES  -10" 

RHOCS = 00 
RHOS  =1.2x10 

Namelist  /SHLDAT/  IB  field  layer  data) 
LCELLS  - 2 

MATERS  - 2.02 

POWS  - 10 

RN  - 15.0 
GAMMAS  = 2 
PRES  = 106 

VN  = 0.0 
RHOCS  = 0 0 
RHOS  =1.2x10 

Namelist  /SHLDAT/  (liquid  (water)  liner) 
LCELLS  - 15 

MATERS  - 3.03 

POWS-  2.0 

RHOS  — 1. 000000000 1 

RN  = 30.0 
GAMMAS  = 0.5 
PRES  = 10h 

VN  = 0.0 
RHOCS  = 10 

Namelist  /SHLDAT/  (piston  layer) 

LCELLS  - 5 

MATERS  - 4 04 

POWS  = 1 5 

RHOS  - 7 80000000078 

RN  - 35.0 
GAMMAS  = 0.5 
PRES  = 106 

VN  = 0.0 
RHOCS  = 7.8 

Namelist  /SHLDAT/  (driver  gas  plenun) 
LCELLS  - 5 

MATERS  - 5.05 

POWS  - 10 

RN  = 40.0 
GAMMAS  = 1.4 
PRES  =10* 

VN  = 0 0 
RHOCS  = 0.0 
RHOS  = 0.12 

Several  outputs  from  the  complete  calculation  are  included  in  Appendix  E.  A number  of 
ad  hoc  modifications  of  the  test  program  were  also  included  for  this  particular  calculation.  In 
the  timestep  control  portion  of  the  program,  the  initial  limestep  is  taken  as  the  geometric  mean 
of  the  minimum  and  maximum  values  specified.  In  a problem  where  geometric  convergence  is 
strong,  as  in  standard  test  #3,  the  spread  between  8rmm  and  8rmav  has  to  be  large  Thus  the  ini- 
tial timestep  is  unspecified  when  all  the  velocities  are  zero. 

This  type  of  calculation  starts  with  an  implosion.  The  liner  accelerates  radially  inward 
until  the  plasma-magnetic  field  core  becomes  highly  compressed  When  the  pressure  interior  to 
the  liner  is  high  enough,  the  implosion  comes  to  rest  and  is  then  converted  to  an  accelerating 
explosion.  At  "turn  around"  the  liner  comes  to  rest.  At  this  time,  the  timestep  computed  by 
DTFLOW  is  large  because  the  fluid  is  moving  slowly  even  (hough  the  accelerations  are  max- 
imum. To  avoid  the  problem  of  taking  too  large  a timestep.  users  might  be  well  advised  to 
include  in  the  timestep  control  portion  of  their  codes,  an  estimate  of  the  change  in  velocits 
expected  during  the  step  as  well  as  the  actual  velocity  at  the  beginning  of  the  step  In  the  test 
program  here  this  problem  has  been  tackled  in  two  ways.  Both  the  old  velocity  and  the  current 
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LINUS  SIMULATION  INITIAL  CONDITIONS 


INITIAL  VELOCITIES  ZERO 


f igure  10  — A schematic  of  the  LINUS  simulation  initial  conditions  A five  laver  calculation  involving 
liquids,  gases,  and  plasma  Parameters  of  the  calculation  are  close  to  those  used  in  recent  LINUS  I experi- 
ments at  NRl  with  the  exception  of  replacing  the  plasma  and  magnetic  field  with  a ft  - 14  gas 


velocity  are  used  to  get  the  smallest  estimate  of  8/  available.  Then  the  estimated  5 r is  limited  to 
be  no  more  than  1.1  limes  the  previous  limestep  This  last  operation  ensures  that  the  step 
doesn’t  automatically  become  large  when  the  liner  comes  to  rest. 

The  test  program  also  has  a section  of  code  at  the  beginning  of  the  diagnostics  portion 
uhich  reduces  the  print  out  interval  near  liner  turnaround.  Since  this  occurs  after  cycle  650, 
the  special  tests  have  been  left  in  the  main  program  The  other  tests  performed  for  problem 
#1  and  problem  #2  generally  require  less  than  a couple  ol  hundred  steps  so  this  problem 
dependant  code  is  never  invoked  During  the  deceleration  and  turnaround  period  the  explicit- 
ness parameters  *,  and  t,  are  also  set  to  zero,  an  ad  hoc  fix  for  the  nonlinear  instability  and 
non-convergence  problems  which  arise  when  a centered,  non-dissipalive  algorithm  is  subject  to 
fast  transients.  The  limitation  here  seems  10  be  the  fact  that  the  solid  equation  ot  stale  cannot 
permit  p < p,,  a situation  which  occurs  behind  the  water  hammer  unless  strongly  damped 

higure  11  displays  the  radial  position  of  the  inner  edge  of  the  linear  as  a function  ol  time 
near  turn  around  Some  asymmetry  seems  evident -possibly  correlated  with  the  weak  compies 
sions  of  the  liner  material  noted  near  turnaround  in  the  appendix. 
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Figure  !!  — Trajectory  r vs  / of  the  inner  edge  of  the  liner  at  turnaround.  A very  slight  compression  of  the 
liner  material  is  expected  for  these  parameters  with  a corresponding  modification  of  the  trajectory  near  tur- 
naround 


VII.  DISCUSSION 

ADINC  was  originally  motivated  by  four  related  but  distinct  computational  problems. 
The  first  of  these  was  the  need  to  extend  our  detailed  reacting  shock  model17  to  treat  flame  pro- 
pagation and  similar  reacting  flow  problems  where  chemistry  and  fluid  dynamics  interact  on 
timescales  appreciably  slower  than  the  sonic  transit  time.18  The  self  consistent  energy  release  of 
many  reacting  flows  occurs  in  extremely  narrow  reaction  zones.  Therefore  Eulerian  algorithms 
are  dangerous  because  their  unavoidable  numerical  diffusion  u'lf’  usual  swamps  real  molecular 
diffusion  unless  the  grid  resolution  is  unacceptable  fine.  Thus  the  Lagrangian  aspect  of  the  cal- 
culation  is  also  important.  We  could  have  used  the  ID  analog  of  our  2D  slow  flow  algorithm5 
for  many  of  the  problems  but  felt  that  the  ability  to  follow  sound  waves  as  well  would  be 
important  in  studies  of  the  deflagration  to  detonation  transition. 

The  second  problem  which  motivated  ADINC  was  the  potential  importance  of  acoustic 
waves  (weak  shocks)  and  modest  compressibility  in  the  physics  of  imploding  liner  systems 
(LINUS'1'8)  as  discussed  above.  Sonic  transit  times  in  the  liquid  metal  liner  are  up  to  two  ord- 
ers of  magnitude  shorter  than  the  implosion  lime  thus  implicit  integration  is  again  important. 
Furthermore,  the  large  density  discontinuities,  coupled  with  the  diverse  materials  and  discon- 
tinuous properties  in  even  the  simplest  LINUS  system  clearly  demand  a Lagrangain  treatment 
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of  at  least  the  material  interfaces. 

The  third  problem  which  led  to  ADINC  was  the  desire  to  develop  a fully  compressible 
version  of  our  two-dimensional  hydrodynamics  code  SPLISH 19  20  which  would  allow  discontinu- 
ous material  interfaces  and  reasonable  approximations  to  real  equations  of  state.  SPLISH  is 
Lagrangian  and  uses  a reconnectable  grid  of  triangles  to  permit  long  time  integration  of  strongly 
sheared  flows.  Because  of  the  complexity  of  the  SPLISH  variable  geometry,  it  was  clearly  going 
to  be  profitable  to  perform  tests  of  the  proposed  compressible  algorithms  in  one  dimension  first. 
ADINC  fills  this  role  as  well. 


j 
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The  fourth  problem  which  motivated  ADINC  is  the  detailed  simulation  of  a laser-driven 
ablation.  The  detailed  analysis  of  the  structure  of  this  ablation21  predicts  a thickness  for  the 
ablation  layer  of  less  than  0 1 n but  the  sound  speed  in  this  region  can  exceed  I07cm/sec. 
Thus  picosecond  timesteps  would  be  required  for  explicit  courant  stability  — a strong  induce- 
ment to  look  for  accurate  implicit  algorithms  since  10  nsec  and  longer  experiments  have  to  be 
simulated.  Here,  as  in  the  LINUS  problem,  a single  calculation  must  span  the  range  of  tensi- 
ties from  solid  or  greater  to  a tenuous  gas.  Furthermore,  the  vast  span  of  zone  sizes  required 
for  different  regions  of  the  flow  clearly  requires  special  treatment  to  retain  accuracy. 

The  three  test  problems  presented  above  were  designed  to  give  the  prospective  user  of 
ADINC  a useful  background  of  experience  with  the  program  in  different  regimes  and  witn 
different  problems.  Various  types  of  errors  were  tested  and  discussed  but  clearly  a lot  of  work  is 
required  on  shocks  and  other  strong  nonlinear  sonic  phenomena.  In  particular,  until  an  adap- 
tive rezoning  algorithm  is  included,  even  the  flexible  initial  zoning  allowed  by  varying  the 
parameter  POWS  in  Namelist  /SHLDAT/  is  not  adequate  to  deal  accurately  with  a number  of 
problems.  Auxiliary  reports  in  this  series  will  present  and  test  additions  to  the  ADINC  package 
for  this  rezoning,  for  more  realistic  equations  of  state,  and  for  propagating  discontinuities  and 
shocks. 

ADINC  is  an  evolving  package  and  its  users  are  part  of  the  evolution  process.  Your  com- 
ments and  suggestions  for  improving  or  simplifying  the  calculation  will  certainly  be  taken  into 
account  in  future  editions  of  the  package.  Comments  on  omissions  and  inaccuracies  in  this 
documentation  will  also  be  gratefully  received  and  incorporated. 
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Appendix  A 

THE  ADINC  PACKAGE  SUBROUTINES 


r 

l 


r 


1 

2 

3 C 
H C 

5 C 

6 C 

7 C 
9 C 
9 C 

10  C 

11  C 

12  C 

13  C 
la  C 

15  C 

16  C 

17  C 
16  C 

19  C 

20  C 

2 1 C 

22  C 

23  C 
2a  C 

25  C 

26  C 

27  C 
29  C 

29  C 

30  C 

31  C 

32  C 

33  C 
3a  C 
35  C 
3b  C 
37  C 
39  C 
39  C 
ao  C 
«1  C 
«2  C 
u 3 C 
aa  C 
«5  C 

«b  c 

«7  C 

ae  c 

ap  c 
5o  c 


SUBROUTINE  ADInC  (raD,  VEL » RHo,  PRE,  N,  OTIN,  CYCLE, 
t RRNtw,  renew,  vRNEw*  vl.NEW) 


aoinc  has  been  constructed- as  a utility  package  to  advance  thc 
FOUR  HYORODYNAHlc  variables.. 

RAD(I)  * POSITION  (RADIUS)  OF  THE  I-TH  CELL  INTERFACE  C C J 

VE L c 1 3 8 VELOCITY  OF  THE  1-th  CELL  INTERFACE  Ccm/SEC) 

RHO(I)  » DENSITY  in  CELL  I between  INTERFACES  1,1+1  (GN/CC) 

PREin  8 PRESSURE  IN  THE  I-TH  COMPUTATIONAL  CELL  (ERG/CC) 

lagpangian  fluid  dynamics  equations  are  s*lvec  including  a flex- 
ible equation  of  state  which  can  vary  from  cell  to  cell  ifi  the 
DISCRETIZED  REPRESENTATION  OF  the  FLUID.  The  EQUATIONS  SOLVEr  ARE 


D (RAD)  D (VEL ) - 1 

......  ■ VEL,  ......  s .....  c,raD  (PRE), 

DT  DT  RhC 


AND  THE  EQUATION  Op  STATE  .. 

( PRF  ) *«  1/GAMMAC 
RHO  « RHOC  ♦ (-....) 

( ENTC ) , 

all  NON-IDEAL  EFFECTS  which  might  re  INCLUDEC  in  Af  ATPC  CALCU- 
LATION have  TO  5t  INCLUDED  SEPARATELY  EITh^p  pv  PHENOMENOLOGICALLY 
IMBEDDING  A SIMPLE  MODEL  IN  The  CALCULATION  OR  py  TIMESTEp 
SPLITTING. 


EACH  FULLY  LAGPANGIAN  CELL  HAS  SEVERAL  QUANTITIES  THfcT  ARE 
CONSERVED  MOVING  with  THE  fluid  AS  LONG  AS  diffusive  AND  OTHER 
NOn-IOEAL  EFFECTS  AND  SOURCE  TERMS  APE  NOT  INCLUDED  In  ThE  CALCU- 
LATION. THE  EQUATION  OF  STATE  IN  EACH  FLUID  CFLL  maY  DIFFER. 

THE  QUANTITIES  INVOLVED  In  THE  EQUATION  OF  STATE  a RE  INITIALIZED 
by  the  two  ENTRIES  SETMat  and  SETEOS  AND  The  equation  «f  statf  is 
evaluated  by  calling  useeos.  thf  egua t i of  of  statf  quantities 
ape  COMMUNICATLD  Throughout  the  aoinc  Package  in  common  BLOC* 
/adicom/.  these  "Constants"  vary  from  cfll  to  cell  according  t* 
the  initial  COf'DlTiONS.  FOLLOWING  are  ThE  DEFINITIONS  of  These 
QUANTITIES. . 


MATEKCCI) 


MASSC(I) 
gammac  ci ) 
ENTC(I) 
P-HOC  (1) 


• CELL  IDENTIFIER  o L.mm  WHERE  n < L < in  IS  The 

LAYER  MIMpFP  Af'C  C < mm  < loo  IS 
THF  material  identifilr. 

s CELL  MASS  e RHO(I)*LAM(I)  - HELD  CONSTANT  IN  ADINC 
e cell  ADIABATIC  r.AS  CONSTANT  - HF L D F I XFD  IN  AT  I f.C 
E cell  ENTROPY  - CONSTANT  DURING  ADUC  HYC-ROD  YN  4m  i c 0 

* DL,!SITY  CONSTANT  IN  THE  EQUATION  or  STATE  CGM/CC) 
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■ 


V- 


V 


51  C nctlis  * M'MFiEP  OF  CELLS  OF  FLUID  I N The  CALCULATION 

52  C 

5 3 C V4FIABLE  C- A m P a and  entropy  ape  USED  in  EACH  LaGSangTaN  cell. 

54  C Af  IMPLICIT  FPESSL'PF  ITFRaTION  ENSURES  LINEAP  STABILITY  BUT  high 

55  C FREQUENCY  PHEN0REnA  are  INACCURATELY  INTEGRATED  WHE  N THE  TimeSTEPS 

5b  C ACE  CHOSEN  Tr  E’F  APPRECIABLY  LONGER  THAN  The  COURanT  TImeSTEP.  The 

5y  C NONLINEAR  TERMS  APE  ITERATED  WITH  A QUADRAT ICALLY  CONVERGENT  ALGO- 

58  C RlTHfi.  PEFERENCEi  URL  MEMORANDUM  REPORT  «????.  1979. 

5o  C 

60  C PROBLEMS  IN  ORE  OF  POUR  GEOMETRIES  CAN  PE  SET  Up  FOR  ACINC  BY 

61  C CHANGING  the  integer  alpha  in  the  Call  to  SETGEO.'.. 

6?  C 

65  c alpha  = i cartesian  coordinates. 

64  C alpha  s 2 CYLINDRICAL  COORDINATES. 

65  C ALPHA  b 3 SPhep  ICAL  COORDINATES. 

6b  C alpha  B 4 P1*Ep  SERIES  COORDINATES. 

07  C 

66  c AOpC  USES  The  UTILITY  USFGEe  to  DETERMINE  the  INSTANTANEOUS  GPID 

69  c quantities,  call  usegeo  calculates  the  interface  areas,  the  cell 

To  C VOLUMES,  AND  THE  CELL  CENTER  POSITIONS.  ADINC  DOES  NOT  UPDaTE  ALL 

7 l C the  GEOMETRY  QUANTITIES  AUTOMATICALLY  ON  EXIT.  THUS  USEGE®  must 

Ta  C aLS"  BF  called  EXTERNALLY'.  IN  GENERAL  ADINC  DOES  NOT  HAVE  ACCESS 

73  C TO  the  USER-DEFINED  ARRAYS  FOR  AREA,  RADC,  AND  LAMC. 

74  C 

75  C THE  BOUNDARY  CONDITIONS  TREATED  BY  aDInc  APE  QUITE  GENERAL. 

76  C THE  POSITIONS  AND  VELOCITIES  *F  THE  REGION  BOUNDING  INTERFACES 

77  C PAO(l)  AND  RAD(N+1)  CAN  BE  EXTERNALLY  DETERMINED  FUNCTIONS  OF  TIME 

78  C AND  OTHER  PHYSICAL  VARIABLES  DURING  THE  CALCULATION.  THE  BOUNDARY 

79  C CONDITIONS  APE  CChNUNICATED  TO  THE  ADINC  PACKAGE  VIA  THE  POUR 

BO  C ARGUMENTS... 

01  C 

82  C RRnEW  « RIGHT  BOUNDARY  POSITION  RAD(Ntl)  AT  END  OF  TIMESTEP 

8 3 C RLNEW  * LtpT  BOUNDARY  POSITION  RaDCU  AT  END  OF  TIM£StEP 

84  C VfiNEW  * right  BOUNDARY  VELOCITY  VEL(Ntl)  AT  END  OF  TIMESTEP 

85  C VLnEW  s LEFT  BOUNDARY  VELOCITY  VEL  C 1 J AT  END  OF  TIMESTEP 

86  C 

87  C SEVFPAL  AUXILIARY  VARIABLES  ARE  USED  By  ADINC  ITSELF  OR  THE 

88  C ADINC  ROUTINES  WHIch  SHOULD  ALSO  BE  EXPLAINED  TO  THE  USER. 

89  C 


9o  C 
°1  C 

’?  c 

9j  C 
9a  C 
R5  C 
R 6 C 
97  C 
96  C 
99  C 

100  C 

101  C 

102  C 
K'3  C 
104  C 
1 ('5  C 
1 O0 

1 n 7 
t OB 
1 09 
1 1 0 
1 1 1 
1 12 


N s NCELLS  ■ NUMBER  OF  FLUIO  CELLS  IN  THE  ADINC  INTEGRATION 
DTIN  B T I ME  INTERVAL  FOR  THE  ADINC  INTEGRATION  WHlCN  M*v 

SURCYCLE  UP  TO  100  TIMES  INTERNALLY  IF  NEEDED  FOR 
ACCURACY  OR  STABILITY. 

CVCLE  * T IMESTEP  NUMBER  USED  BY  ADINC  FOR  IDENTIFICATION 

EPSRO  * EXPLICITNESS  PARAMETER  FOR  THE  POSITION  INTEGRATION 

EPSV0  * EXpLICITNESS  PARAMETER  FOR  THE  VELOCITY  INTEGRATION 

N C a m P B NUMBER  OF  CYCLES  AT  THE  BEGINNING  OF  A CALCULATION 

IN  WHICH  ADDITIONAL  DAMPING/SMOOTHING  IS  APPLIED. 
EPSR  ■ EXPLICITNESS  PARAMETER  FOP  RAD  LAST  USED  BY  ADINC 

Epsv  S EXPLICITNESS  PARAMETER  FOP  VEL  LAST  USED  BY  ADINC 

THE  PARAMETER  NPT,  hepE  202,  must  EE  At  LEAST  TWO  LARGER  ThA‘ 
THL  NUMBER  OF  FINITE  DIFFERENCE  CELLS  BEING  INTEGRATED  BY  ADINC. 


FArAMETFR  MPT  = p02,  MPT  s 2*NPT 

INTEGFR  cycle 

REAL*8  RAD(NPT),  VEL(NPT),  RHO(NPT),  PRE(NPT) 

PtAL*8  RliOOOlPT},  RHON(NPT),  pReO(nPT),  PREfUNPT) 

RE  Al  *8  LAMOCNPT),  LAMN(NPT),  RADCO(NPT),  RaDCN(NPT) 

RtAL*8  PIN(NPT),  VIN(NPT),  DLAM(NPT),  LAMEOS(NPT) 

PEaL*8  RBARI(NPT),  PPaR(NPT),  RDRl(NPT),  DTORHO(NPT) 
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1 


113 

PEAL  *8 

AA(')PT), 

PB(MPT), 

CC (MPT) , 

no (mpT) 

1 U 

PC  AL  *8 

A I ( N°T  ) , 

BI(NPT), 

C (NPT), 

OP (NPT) 

1 1 5 

PFAL*8 

DM (NPT)  , 

SCA(MPT), 

SCB(MPT), 

DLAMpP(NpT) 

116 

PE  AL  *8 

AREaO(NPT), 

aREAH(nPT), 

AREAN (NPT) 

117 

REaL*8 

errlIm, 

ERRmaX  , 

ERROR, 

amnmaS 

118 

Pf AL*8 

PREMAX, 

PREMIN, 

RPREMX , 

DT 

11R 

PEAL*8 

V PNE  w , 

VLNEW, 

RPNEW, 

RLNEW 

120 

Rr  al*R 

vpold, 

VLOLD, 

RR  OLD , 

rlOld 

121 

peal*« 

DVPOLC, 

PVLOLD, 

dprold. 

drlold 

122 

REAL*8 

EPSP, 

OMEPSR, 

EPSV, 

OMEPSV 

123 

PE  AL  *8 

EPSP a, 

EpS  VO , 

OMERDT, 

A M X M A S 

124 

REAL  *8 

dtin, 

DTVAL 

125 

C 

12b 

C 

DECLARATIONS 

FOP  COMMON 

BLOCK  /aDICOM/  APPEAR 

THROUGHOUT  THE 

127 

C 

PACKAGE  ROUTINES. 

128 

RE  AL  *8 

MASSC(NPT) , 

ENTC(NPT), 

RHOC (NPT) 

9 

HATERCCNPT) 

129 

PEAL  *8 

gammac (NPT) 

130 

COMMON 

/ADICOH/  HATERC,  MASSC, 

GaMmaC, 

ENTC,  RhOC,  NC 

131 

c 

132 

DATA 

ndamp  / X 0 / , 

epsr,  EPSV 

/O.ODO, 

0 

.oro/ 

133 

DATA 

nc*ll,  niter,  ntihe  /o, 

0,  0/ 

13a 

DATA 

EPRLlM,  EPSRO, EPS VO, ITEMAX/ 1 .OD-9 

$ 

O.aSDO,  0 . "5 

AO  I NC  FORMATS  FOR  DIAGNOSTIC  AND  error  PRINTS. 

format  (U'ADINc  TIMESTEP  PROBLEM  at  CYCLE',  15,  ' DT  ■ 
i 1PD>2.«,  ' ANO  ITEMAX  ■ ',  12,  ' EPRLIN  ■ ',  012.",  /, 

> 1CX,  ' E°RMAX  AT',  I",  ' 8 ',  012.",  ' AND',  13, 

5 ' CELLS  NOT  CONVERGED.') 

FORMAT  C'OADINc  TIMESTEP  PROBLEM  AT  CYCLE',  15,  ' DTVAL  8 ', 
[ 1 PD  1 2. " , ' BUT  DT I N 8 ',  D12«",  ' TIMESTEP  SUBCTCLED.  ' ) 

FORMAT  CfADINc  TIMESTEP  PROBLEM  AT  CYCLE',  15,  ' DTVAL  ■ ', 


1PD12.",  ’ BUT  DT IN  s ',  D12.",  ' CALCULATION  STOPPED.') 
FORMAT  CoaDINc  INPUT  PROBLEM  AT  CYCLE',  15, 

' THE  SYSTEM  SIZE  N I",  ' OUT  OF  RANGE.  NPT  8',  la, 

' CALCULATION  STOPPED. ' ) 

Format  C'OADINc  INPUT  PROBLEM  AT  CYCLE',  15, 

' THE  LEFT  BOUNDARY  HAS  CROSSED  The  RIGHT  .',  1P2D12.4 
' CALCULATION  STOPPED.  ' ) 

FORMAT  ( ' OAR  IN'C  input  proplfm  AT  CYCLE',  15, 


1P2D12.", 


1 ' THE  CELL  SIZE  ^AS  1PD12.",'  at  CELL  ',  I", 

5 ' calculation  STOPPED.  ' ) 

FORMAT  C'PADI'C  INPUT  PROBLEM  AT  CYCLE',  15, 

! ' THE  DENSITY  min  was  ',  1 PD  1 2 . a , ' AT  CELL  ',  la, 

’ ' CALCULATION  STOPPFC.') 

FORMAT  C'OADINc  INPUT  PROBLEM  AT  CYCLE',  15, 

I I the  PRESSURE  him  waS  ',  IP(:12.",'  AT  CELL  ',  I", 

? ' CALCULATION  STOPPFD.  ' ) 

FORMAT  ( ' ('  A D I N r FPEGUFNCY  COUNTERS  (SINCE  LAST  CHECK)  AT  ', 
l 'CYCLE',  15,  /,  1 0 X , 'NO.  CALLS  =',  15,  ' NO.  TlMf, 

? 'STEPS  st,  15,  ' TOTAL  NO.  ITERATIONS  »',  15,  /) 

CHECK  THE  INPUT  TO  ADIHC  FOR  REASONABLENESS. 
fjPTO  s NPT 

IF  C N.LE.2  .OR.  N.GT.NPT-2)  WRITE  (6,  1U0U)  CYCLF,  N,  NPTO 
IF  ( H.LE.2  .OR.  M.GT.NPT-2)  STOP 

IF  (RRNEW.LE.RLNEW)  write  CO,  1005)  CYCLE,  BLNE*.  PRNEw 

if  (PPNEW.LE.pLNEW)  stop 

DO  21)  I ! 1|  N 

OMCI  + I ) s PADClM  ) - BAD  C I ) 

CALI  NAXMIN  (Ur(2),  N,  premax,  IMAX,  PREMIN,  IMIN) 


1 


175 

IF 

( PPE m I N , Lf . O.ODO) 

WRITE  (6,  1006) 

CYCLF, 

PREMIN, 

I M I N 

176 

IF 

(PPEM1N  .lL.  O.ODO) 

STOP 

177 

CALL  MAXHIN  ( rH0 ( 2 ) , N, 

premax,  imAx, 

PREMIN, 

ImIN) 

178 

IF 

(PRLMIH  .lf.  O.ODO) 

WRITE  (6,  1007) 

CYCLE, 

PREMIN, 

I M I N 

17o 

IF 

(PREMIN  .lC.  O.ODO) 

STOP 

1 By 

CALL  m A X m I N ( pRF ( 2 ) , N, 

PREMAX,  IMAX, 

PREMIN, 

I M I N ) 

18j 

IF 

(PREMIN  .LF.  O.ODO) 

wPITF  (6,  1008) 

CYCLE, 

PREMIN, 

1 M I N 

182 

IF 

(PREMIN  ,lE.  O.ODO) 

STOP 

'83 

tea 

tes 

ISt, 
1*7 
1 68 
189 
1 9o 
I'M 

192 

1 Q5 

194 

195 
19b 
197 
19? 

199 

200 
201 
202 

203 

204 

205 
2ob 

207 

208 

209 

210 
2U 
2 1 2 

213 

214 

215 

2 l to 

217 

218 

219 
2?0 
221 
222 

223 

224 

225 

220 


ESTABLISH  INTEGRATION  AND  CONTROL  CONSTANTS  FOR  this  CYCLE. 

M « N ♦ 1 

CALL  maXMIN  CmaSSC(2),  N,  amxmaS,  IMAX,  AMNMa$,  Imi*) 
EPRPfi  s PSQRT (AMXMAS/AMNMAS)*ERRLIM 
RAOCN’ ) 


R9PLC  < 
RL OLD  i 
V»OLD  i 
VLOLD  i 
NCALL  ' 

NCALIH 

EPSR  ■ 
EPSV  ■ 
OHCPSR 
Pmepsv 


pao(i  ) 

VEL(N' ) 
VELCl) 
NCALL  t 


1 


• mjno  (NCALL » NDAHP) 
DELOAT(MCALIM)*EPSPO/DFLOaT(NDAHP) 
DFLCaTCNCALIH)*F.PSVO/CEL0aT(nDAmP) 

• 1 .01)0  - EPSR 
> 1 , fiDO  - EPSV 


30 


CHECK  THE  TIMESTEP  fop  SUPCYCLING  and  NOTE  any  PROBLEMS. 
call  PTEl"W  ( 9 AD # VEL,  DTVAL,  N). 

NSTEP  > 1 
DT  c DTIN 

IF  (DTVAL  ,GT , DTIN)  GO  TO  40 
IF  (DTVAL  . L T , o.01Do*DTIN)  GO  TO  30 
WRITE  (6,  1C02)  CYCLE.  DTVAL,  DTIN 
NSTEP  s (DTIN  + DTVALl/DTVAL 
DT  « DTIN/EL0*T (NSTEP) 

GO  TO  40 

WRITE  (6,  1003)  CYCLE,  DTVAL,  DTIN 
STOP 

INITIALIZE  variables  FOP  the  SUBCYCLING  and  ITERATIONS. 


40 


50 


55 


DVROLP  ■ 

C VL OLD  * 
DPROLD  • 
ORLOLD  * 
DO  50  I 
R IN  ( I ) ■ 

V I N ( I ) « 
RIN(1)  « 

R In (ii i ) * 

V I r.  ( 1 ) * 
vincnd  ■ 
DO  55  1 » 
PREOCI)  * 
P R E N ( I ) 


(VRNEw 
( V L N E W 
(RRNEW 
(RLNEW 
! 2 , N 
R AD ( I ) 

VELCl) 

PLOLD  ♦ ORLPLD 


V»OLD)/FLOAT(NSTEP) 
VL'Tld)/FLPAT  ( N S t E p ) 
RROLD)/FLCAT(NSTEP) 
RLOlC)/FLPAT(NSTep) 

♦ DT  *vEL ( I ) 


RPOLC 
VLOLD  ♦ 


CPPOLO 

DVLOLD 


V R C L f.  ♦ DVROLD 
2, 

FRE  Cl) 

PRE(I) 


227 

CALL 

USE  GE  0 

C p AD , 

aRCaP, 

RADCP, 

LAHO,  N) 

■>28 

CALL 

I’SEEPS 

(WHOP, 

ppfP, 

lameps, 

DLAMOP, 

N) 

2 9 

CALL 

USEGEP 

(»IN, 

abean. 

RADCN, 

LAMN,  H) 

• 

230 

call 

USEEPS 

(RHOf., 

PREN, 

LAMEPS, 

DLAMOP, 

N) 

23  | 

232 

233 

234 

235 

23b 


RPRfl'X  ■ I.ODO/PREMAX 
DO  60  I ■ 2,  N) 

6C  PRECI)  • O.ODO 

PERFORM  TIMESTEP  SUBCYCLING 
CO  550  ISTEP  ■ 1 , NSTEP 
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w 


NT  I Ht  » NtImE  ♦ 1 


237 

236  C 

23”  c PFfjFSRM  TH£  I T£Ra  T J pM  F"P  NEW  VALUES  RHON,  L A MN  , PR£N,  PJN  AND  VIM 

2«0  DO  500  ITER  = 1,  ITEHAX 

2«1  N1TFR  s NITER  ♦ 1 

2<*2  C 

2«3  C CALCULATE  QUANTITIES  USED  IN  TRIDIAGONAL  EXPRESSIONS, 

2<44  C'6  200  I * 1,  Nl 

205  AREA H (I)  s 0 . 5D0» C APE  A 0 C U ♦ AREAN{I)) 

2<*b  2°0  HP API (I ) = <*,5qo*  (RAD  tl ) ♦ R I M ( j ) ) 

297  U o 205  I • 2,  N 1 

2«e  DLAMtn  = LAMrestn  - launch 

2«9  205  RBARCI)  S 0.50C*(RAPC8(n  ♦ RaDCN(I)) 

250  08  210  I * 2,  f. 

25 1 RDPIC!)  * I H0N(I)*(R6aPI  (I)  * PBAR(I))  ♦ PH8N ( I ♦ i ) * ( PB aR (I ♦ 1) 

252  1 - F, B A R I C I ) ) 

255  ?K,  DTCRhOCI)  = PT/POPKI) 

25u  C 

255  C CALCULATE  EXPRESSIONS  USED  'N  THE  TRIDIaGONaL  C 8EEF  I C I ENTS  . 

256  twEPPT  * -ORLPSR*DT 

257  D 8 250  I = 2,  \ 

256  AI(I)  = VELCI)  - DT0Rh8(I)*EPSV*(PRe8(I*1)  - PHEOCI)) 

25”  25o  Bid)  5 DTORM”  d ) *ORFPSV 

2 bo  08  260  1 e 2,  M 

261  Ppd)  e 8f1ERrT*ARFAH(i) 

2b2  DM(I)  s 0NERDT*AREAH(I-1 ) 

265  CCI)  a PREN(I)*DLAMDP(l)  ♦ DPd)*VINCI)  - DM(I)*VIN(I-1) 

26R  260  C(I)  = C(I)  -ULAM(I) 

265  C 

266  C CALCULATE  THE  TPIPIAGCNAL  COEFFICIENTS  WITH  BOUNDARY  CONDITIONS, 

267  DO  500  I * 2,  M 

266  LDCI)  s CCI) 

26”  5(iy  6«d)  = DL  AHdP(J) 

270  DO  505  1*2,' 

271  CCCI)  » -C’PtI)*Bim 

272  AA(I+1  ) * -DN(l  + l)*PId) 

275  bf*Cl)  = E-  6 C I ) - CCCI) 

27a  505  DD  C I ) » DC ( I ) - DP (I) *A 1(1) 

275  CO  510  I * 3 i HI 

270  tiP  c i ) = eeci)  - A A ( I ) 

277  51  0 DD  C I ) = DC'  C I ) ♦ DH  ( I ) * A I C I- 1 ) 

2Tp  t r.  c 2 ) s DO  ( 2 ) ♦ DH  C2)  * VLNEW 

27q  r.D(Nl)  = DDCM)  - DP  (N  1 ) «VRnEh 

2s''  A A C 2 ) = 0 • 0 D o 

28 i CCCN1)  * O.ODO 

2P2  C 

285  C SOLVE  THE  TBI DIAGONAL  SYSTEH  AND  CONSTRUCT  THE  N£ W VALUES  AT  TIHe 

2«a  C t , CT  fob  ThC  CURRENT  ITERATION.  TPIDDV  FAILS  WHFN  n < j5  ThUS 

285  C TRIDOS,  THE  SCaLas  VERSION,  IS  USFD  IN  THIS  PEGIHE  INSTEaD, 

28b  IF  CN  ,LT.  15)  CALL  TRIODS  (N,  A A C2 ) * BE?C2) » CCC2),  CD C2)» 

2e7  1 PPLNC2),  Sc  A C 2 ) , SCO  C2)  ) 

2sg  IF  CN  ,&E.  15)  CALL  TRIDDV  (N,  A A C 2 ) * BBC2),  CC  (2) , DD£2)» 

28R  1 p R L N ( 2 ) , SC  A ( 2 ) , SCB ( 2 ) ) 

290  C 

291  c IN  SOHf  CIRCUMSTANCES  IT  may  be  APPROPRIATE  T*  LET  thf  PRESSURE  GO 

292  c NEGATIVE,  IS  SUCH  CASES  The  FOLLOWING  LOOP  HUST  DF  °Emovf T OR 

295  C EREMIN  MUST  BE  ALLOWED  TO  GO  NEGATIVE. 


29a 

DO  340 

I 

= 2,  M 

295 

PRENCI) 

s L H A X 1 

(0 .0 IoOaPREMIN,  pre* 

CD) 

c90 

L’O  350 

I 

= 2,  N 

297 

V I N ( I) 

s 

AICI)  - 

eid)«CPREN(I*l)  - 

PRENCI) ) 

29b 

3*5  0 

TN(I) 

= 

Fi At  (I) 

♦ LT*CEPS»*VLLCI)  ♦ 

omepsr*v  In  CD) 
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call  I'SEGrO  CbIN,  aPL  an 1 R ADC  n , lawn,  N) 
call  USEEOS  ( R H o N , PRf'J,  LAIC'S,  PLaNDP,  N) 

CHECF  ON  mtlE  THC P THE  ITFPaTI"N  HAS  CONVERGED. 

po  pqo  i s a,  1. 1 

SCA(I)  = Pffrtn  - PRF(I) 

SCACI)  * L'aBsCsCA(I))*EPPFNX 
SCCtl)  * H'M)« 

PRECI)  * PREN(I) 

NOTCCf  * o 
ICLLL  = « 

ERR.MAX  * C.OL'i 
C"  « t! *»  I = ?.  M 

IE  (SCACI)  .LT.  SCT ( i 5 5 GO  T9  405 
N"TCCN  « DPTefti  ♦ 1 
IE  CSCA(I)  .gT.  EPPMAX)  ICELL  s I 
ERRmax  * L’l-Axl  CERP^AX,  SCA(I)) 

CONTINUE 

IF  (NOTCON  .EU.  0)  GO  T 9 505 
CONTINUE 

FR I NT  OUT  Ip  WL  mAVE  NOT  cO'VFPGEP. 

IE  (KOOCCVCLE*  IP)  . FD . 0) 

WRITE  (6,  lOul)  CYCLE,  PT,  ITEf'AX,  EPpOR,  ICELL.  FRRNAx, 
NOTCON 
CONTINUE 

SET  UP  F op  ANOTHER  SUBCVCLE  TIRESTEP. 

IF  ( I STEP  .EG.  NSTEP)  G 9 TP  550 
E9  52"  1=1.  M 
AREAS ( I ) = ATEaNCI) 

VEL(I)  * V I N ( I ) 

RACC1)  = RINCI) 

DO  525  I * 1.  M 

PTNCI)  » EAD(I)  ♦ CT*VtL(!) 

VINCI)  = VEL  CD 
Rif. Cl)  ■ RAO  ( 1 ) ♦ C.  P L * L P 
P I N ( N 1 ) = PAfCM  ) ♦ ORE  OLD 
VIN(I)  « VEL  C 1 ) ♦ DVLOLC 
VIN(M)  * VFL(NI)  + PVROLD 
DO  5 JO  I = 3,  M 
fiHO(I)  S RHONCn 
PRLO(I)  a FPECI) 

RHOOCI)  = RWO^CI) 

CALL  USEGEO  (Pin,  apean,  rapcn,  LAf'N,  N) 
call  U S E E 0 S CrhON,  preN,  LAMEOS,  DLANDP,  n) 

Continue 

cle*n  up  FOR  exit. 

CO  600  I « 1,  HI 
VELtn  * VINCI) 

RAr(i)  » S'  IN  c I ) 

DO  601  I * «?,  Nt 
RHC(I)  * RwOt  en 

RETURN 


E^TRY  APINCO  (NODE,  CYCLE) 


APINCO  PRINTS  CUT  the  NUf'PER  OF  CALLS  TO  AOINC.  THE  OF 


361  C 

362  C 
36  J C 

364  C 

365  C 

366  C 

367  C 
366  C 
366 
37U 

371 

372 

373 

37u  C 

375  C 
37b 

377  C 

376  C 
379  C 

390  C 

391  C 

392  C 
3«3  C 
39a  C 
395  C 
39b  C 
3«7  C 
39P,  C 
399 

3<»0  C 

39  1 

392 

393 

394 

395 

396 

397 
396 
399 


TlMESTfcPS  INCLUDING  SLECVcLING  PERFORMED  BY  ADI^C  DURING  THOSE 
calls,  AND  THE  TnT al  HUHflrR  OF  ITERATIONS  fepforitc  SINCT  THE  last 
call  to  adinco. 

MOoE  » 0 print  out  Aflr  PESET  THE  COUNTERS  IN  A D I NO  • 

node  .ne.  0 beset  all  but  ncall  counter  fob  filtering  operation 

CYCLE  * TjHj-steP  NUHPEP  FOR  IDENTIFICATION  PURPOSES 

HRlTE  (6,  1 0 o 9 j CYCLE,  NCALL,  RTihE,  nitep 

IF  (HOPE  .EQ.  0)  NCALL  s 0 

NT  I HE  s 0 

NITER  s 0 

return 


ENTPV  SETePS  (ERO,  EVP,  MdAmF,  EPPUT,  E VOUT  5 


SFTEPS  PERMITS  The  uSEF  TO  RESET  THE  F.XPLIcITNFSS  PaPaheTERS  'OR 
DIFFERENT  TYFES  OR  stages  OF  fluid  PROllLHS.  the  EXPLICITNESS 
parameters  ape  happed  into  the  fange  0 <*  eps  <«  t. 


EPfl 
E Vo 

M D A H F 

EP3UT 
E V OUT 


s THE  NE R POSITION  EXPLICITNESS  paPAuETEP 
s THE  new  VELOCITY  EXPLICITNESS  PaBametEP 

b the  new  value  of  ndahf,  « or  damping  cycles 
b the  MOST  PFCEMT  P0S I T 1 91.  FXPLirlTNESS  paPamETeP 
b THE  MOST  PF  CENT  VELOCITY  EXPLICITNESS  PaPa^'ME0 


RE  AL *8  ERO , EVO,  EROUT,  EVOUT 


NC AMP  s H A X C (I'D AMP,  1) 

EPSPO  = (MINI  (fro,  l'.ODO) 
EPSVO  s ON INJ  (EVO,  t.GOO) 

E P S R 0 = C m A X 1 (EPSFP,  C.ODO) 
F. P S V 0 s DMAX1  (ERSVO,  o.OOO) 
EROUT  a EFSP 
EVOUT  b EpSV 
RETURN 

END 
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1 

2 

1 

4 

5 

6 

7 

8 
9 

10 

11 

12 

15 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 
2b 

27 

28 

29 

30 

31 

32 

33 

34 

35 
3b 

37 

38 

39 

40 
«1 
«2 
43 

U4 

«5 

«b 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 


1001 


c 

c 

to 


c 

c 

2o 


c 

c 

30 


SUBROUTINE  SeTGIC  (alpha,  GE  ?MC  8 ) 


SETGEO  CONTROLS  THE  GEOMETRIC  ASPECTS  OR  AN  ADINC  INTEGRatIO'  . 
SETGEO  IS  CALLED  6v  THE  USEP  TO  TELL  THE  ADINC  P AC*1  AGE  EXACTl* 
WHICH  OF  TtlF  FOUP  POSSIBLE  GEOMETRIES  he  WISHES  TO  USr. 

ALPHA  s 1 CARTESIAN  COORDINATES. 

ALPHA  ■ 2 CYLINDRICAL  CflOPDINATfS. 

alpha  ■ 3 SPhe R I c Al  COORDINATES. 

ALPHA  a u POW(R  SERIES  COORDINATES. 

GF0MCOC1..S)  « aRHaV  Or  FIVE  COEFFICIENTS  IN  THE  EX°PESSI"N  rOR 
CELL  AREA  AND  INTEGRATED  VOLUME  USED  BELOw. 


SETGEO  INITIALIZES  THE  QUANTITIES  GALPHA,  HALPHA,  G(I),  AND  H C I J 
for  LATER  REPEATEC  USE  IN  ENTRY  USEGEO  BELOW. 

PARAMETER  NPt  a 202 
INTEGER  alpha,  aL.phh 

RE  AL*8  GALPHA,  HALPHA,  G(5),  m(5) 

REAL*8  GE0mC0(5),  R I AM  1 ( NPT ) , CELr(NPT),  TvOL(NPT) 

RE  al  *8  DTH  in (nPT ) , a6SL (NPT) , DPhIN,  Rl 

EQUIVALENCE  (DTMIN(l),  TV0L(1))»  CA9SV(1)»  DELR  M ) ) 

EQUIVALENCE  CRIAMI(I),  TVOL(l)) 


CHECK  the  INPUT  to  SETGEO  AND  INITJalIZF. 

IF  (ALPHA.  L T . 1 .OP.  ALPHA. GT. 4)  WRITE  (6,  1001)  AlPHA,  G£  OMC  0 

IF  (ALPHA.LT. 1 .OR.  ALPHA. GT. 4)  STOP 

FORMAT  C'OSETGeO  INPUT  PROBLEM.  ALPHA  OUT  OF  R a f Gr , 

14,  2X,  1 P5D 12.4) 

PI  « 3. 14)5926535897900 
ALPHH  * ALPHA 

GO  TO  (10,  20i  30,  40),  ALPHH 


HOCE  a 1 RESET  ADINC  FOR  CARTESIAN  COORDINATES. 

GALPHA  * 1.CU0 
HALPHA  * 1.000 
RETURN 


mOcE  a 2 

GALPHA  * 
HALPHA  « 
RETURN 


RESET  ADINC  FOR  CYLINDRICAL  COORDINATES. 
PI 

2.0DC'*PI 


MODE  a 3 
GALPHA 
HALPHA 
RETURN 


RESET  aD INC  POP  SPHERICAL  COORDINATES, 
a 4.cdo»PI/3.OD0 

a 4.0[P*PI 


MOqE  a 4 RESET  ADINC  FOR  POWER  SERIES  COORDINATES. 

DO  41  I a t,  5 
GCI)  a GEOMCCCI) 

MCI)  a 5CI)/DFL0ATCI) 
return 


ENTPY  USEGEO  (RAD#  AREA,  RADC,  LAMC,  N) 


GIVFN  A MONOTONICALLV  INCREASING  SET  OF  CELL  INTERFACE  POSITIONS, 
THE  INTERFACE  areas,  CELL  CENTER  LOCATIONS,  and  cell  volumes  ARE 
CALCULATED  in  a fully  VECYORIZED  manner.  This  geometric  UTILITY 
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r 


f 


65 

C 

IS  i 

USED  BY  ADlNC,  DIAGNOSTICS  ROUTINES  AND  THE  MAlN  PPOGRAM 

m 

64 

C 

WHENEVER  THE  CELL  INTERFACE  CONFIGURATION  is  changed  - TO  UPDATE 

65 

C 

the 

GEOMETRIC  QUANTITIES. 

66 

C 

67 

C 

RADII)  « POSITION  of  THE  I-TH  INTERFACE  Cl  * 1 » N*i) 

CCM) 

66 

c 

aRE*(I)  ■ AREA  IN  the  COMPUTATIONAL  domain  OF  the  I-Th 

CELL 

69 

c 

Interface  ccm**2) 

7 0 

c 

RAOC(I)  •.  POSITION  OF  THE  I-TH  CELL  CENTER  Cl  * 2,  N+l)  (CM) 

7 1 

c 

LAMC(I)  ■ VOLUME  OF  CELL  I BETWEEN  INTERFACES  I,  I ♦ 1 . 

C C M * * 3 ) 

72 

c 

N 

■ NUmBER  OF  INTERIOR  CELLS  IN  THE  SYSTEM 

7 3 

c 

74 

RE  AL*8  RAD(NPT),  AREA (NPT) » RAOCCNPT),  LAMC(NPT) 

75 

c 

76 

c 

77 

c 

CHECK  THE  INPUT  TO  USEGEO  FOR  REASONABLENESS. 

78 

NPTO  ■ NPT 

79 

IF  CN.lE.1  .OR.  N.GT.nPT-2)  WRITE  C6,  1002)  N,  NPTO 

80 

IF  (N.lE.1  .OR.  N.GT.nPT-2)  STOP 

81 

DO  50  I ■ 1,  n 

82 

50 

TVOLCI+1)  = RADCI+1)  - RAO(I) 

83 

CALL  MAXMIN  ( T yOL  ( 2 ) » N,  TV0LC1),  JmAx,  DRMIN,  IMIN) 

84 

if  (DRMIN  ,LE.  0.000)  WRITE  (6,  1003)  DRMIN,  IMIN 

85 

IF  (DRMIN  ,LE.  O.OCO)  STOP 

86 

1 002 

FORMAT  ( 1 OUSE GeO  INPUT  PPOBLEM.  M OUT  OF  RANGE.  14, 

14, 

87 

1 

1 

' calculation  STOPPED. ' ) 

88 

1003 

format  C'OUSEGeO  input  PROBLEM.  CELL  size  NEGATIVE 

89 

\ 

1 

1 PD  I 2 . 4 , ' AT  CELL  14,  ' CALCULATION  STOPPED.') 

90 

c 

91 

NP  a N ♦ 1 

92 

GO  TP  C * 00 , 200  # 300,  4C0),  ALPHH 

’J 

c 

9 a 

too 

DO  101  I a l,  t„p 

95 

toi 

RlAMl(I)  s 1,'>do 

*6 

GO  TO  500 

97 

c 

98 

200 

DO  201  I * 1 , NP 

99 

20  i 

RIAM1CI)  a PAn ( I ) 

100 

GO  TO  500 

101 

c 

102 

300 

DP  301  I * 1,  NP 

103 

301 

RIAMl(I)  a RADCI)«PAD(I) 

104 

c 

105 

c 

FOR 

the  REGULAR  GEOMETRIES  CALCULATE  THE  4PEA  AND  VPLUUE. 

106 

500 

DO  501  I ■ 1,  NP 

107 

AREA ( I ) a HALpHA*PIAMl  (J) 

108 

501 

TVOLCI)  * GALpHA»RIAMl (I)*RADCI) 

109 

GO  TO  600 

no 

c 

in 

c 

for 

The  POWEP  SERIES  (NOZZLE)  COORDINATES  The  USER  UUST  SpFHIFy 

112 

c 

all 

OF  THF  GEOMETRIC  COEFFICIENTS  VIA  ThE  INITIALIZING  AP»ay 

113 

c 

GF 0 Mc 0 IN  ThE  CALL  to  SETGEO.  here  G ( 1 ) S G E 0 M C 0 ( 1 ) » ETC. 

1 t 4 

c 

AREA(P)  * Gl  ♦ G2*P  ♦ G3*R*«2  ♦ g4*R**3  ♦ g5»R»*4 

115 

c 

TVOL(R)  a HUP  ♦ H£*P**2  , h3»R**3  ♦ Hii*R**u  ♦ h5*R**5 

116 

400 

CO  401  I » 1/  NP 

117 

ARfcA(I)  ■ G(5)*RAD(I) 

118 

TVOLCI)  a H,(5)*PAD(I) 

119 

ARE  A ( I ) a PAD(I)*(AREa(1)  ♦ G(4)) 

120 

TVOLCI)  « RAD(I)*CTVOL Cl)  ♦ H(4)) 

121 

AREA(I)  « RAPCI)*(AREa(I)  a G C 3 ) ) 

122 

TVOLCI)  a RAC(I)*CTVOL(I)  ♦ H ( 3 ) ) 

123 

AREACI)  a RA['CI)*(ARE4(I)  ♦ G (2)  ) 

124 

TVOLCI)  « RAJ  CI)*CTV0LCI)  ♦ H(2)) 
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1 


1 


! 


125 

126 
127 
126 
129 
l 30 
1 31 

132 

133 
139 

135 

136 

137 

138 

139 

140 

141 

142 
U3 

1 44 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 
163 
1 64 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 
1 6o 
181 
182 

183 

184 


ABEACi:  s API  A<n  ♦ G(l) 

401  tvol ( i ) * ral (i)*crvoi  m ♦ nun 

c 

c COMPUTE  THE  CELL  VOLUME  AND  CELL  CENTER  LOCATIONS. 

600  CALL  maXMIU  (APEa(I),  HP , K ADC  C 1 J . IMAX,  PRMIN,  jmjkj 
IF  (CP"IN  .LE.  0.000)  wp j Te  ( 6 , 20C 1 ) DRi'IN,  INI' 

IF  (CPMlN  .LL.  0.000)  STOP 

2oO  1 FORMAT  COUStGEO  PRCPlEM,  NEGATIVE  AREA  1P012.4, 

1 ' A 7 CELL  ',  14,  ' CALCULATION  ST*rPEP.') 

CO  601  I * 2,  UP 

LAMCCI)  ■ TVCL(I)  - TvUII-l) 

601  PAOC(I)  « (RAO(I)«APCt (1-1 ) + RArCl-l)*ARFA(I))/ 

1 ( ARE  A C I ) ♦ AREA  C I - 1 )) 

RE TURK 
C 
C 

ENTRY  UTFLOW  (R0u>  vel,  OTVAL,  N) 

c 

c 

C DTfLOW  CALCULATES  A PERMISSIBLE  TIMESTEP  OTVAl  GIVE’  Tl  F SET  or 

C N + l CELL  IN7EPFACES  A NO  THEIR  VELOCITIES. 

C 

c rocci) 

c vr l f i ) 

C DTVAL 

c 

C N 

c 

PE  AL  *8 

c 
c 

c CHECK  THE  INPUT  to  0TFL8W  FOR  REASONABLENESS. 

NPT8  e NPT 

IF  fH.LE.2  ,0P.  N.C-T.nPT-2)  RPITE  (6,  1"04)  N,  NPT* 

IF  (N.LE.2  ,pP.  N.GT.mPT-2)  stop 
Nl  s N ♦ 1 

EPS  b 0 .4900 
op  715  I * 2,  Ml 

715  L ELP  ( I ) » POC’(I)  - ROC(l-l) 

CALL  MAXMIU  (L.'ELRf2),  N,  OELRfl),  ImAx,  r R 1 r l , I*TM 
IF  COPHIU  .LE.  n.CDP)  WRITE  (b,  1005)  DPMIN,  IH' 

IF  (PRMIN  .LI.  r.OCO)  STOP 

1004  FORMAT  I'OCTFLOW  ItiPUT  PEOPLE H.  K OUT  OF  RANGE.  )4,  )4, 

1 1 CALCULATION  STOPPEC.') 

1005  FORMAT  ( 'OOTFLOw  INPUT  PROBLEM.  CELL  SIZE  NEGATIVE 

1 1 PD  1 2 . 4»  ' AT  CELL  ',  14,  ' CALCULATION  STOPPED.') 

C 

c PEqUIRe  VFL*DT  < minimum  of  CELL  WIDTHS  CEL p . 

DO  720  I * 2»  N 

720  DTMIN(I)  b EPSadMINI (rElR(I),  ft  L»  C I ♦ 1 ) ) 

OTHIN(l)  * Cttpt2) 

DTKIN(Nl)  P r£LP(N)  ) 

DO  725  I * 1,  M 
725  APS V ( I ) b [ A [ . S ( v E L ( I ) ) 

DVSAFF  * 1.0D-4? 
no  73o  I « 1,  IJ1 

736  DTMIMI)  « OT"IH(I)/(ABSV(I)  ♦ PVSAFE) 

Call  MAXMIU  (OTMIN,  Nl,  DTMaX,  IMAX,  DTVAL,  I'TO 
RETURN 

END 


* PCSiTior  Of  the  I-TH  cell  INTERFACE  CCu) 

S velocity  of  THE  1-th  CELL  Iutfofacp  (Cy/SEC) 

- THE  ESTI-ATEC  VAl'JE  "F  A PERMISSIBLE  T T M F S T £ P whICh 
PEEVEfTs  I'TFPFACE  CROSSING  assuming  the  "*Tlou  G I v e 4 
b THE  NUMBER  OF  INTERIOR  CELLS  IN  tmf  CUPPFnT  SyStF'’ 

ROD(N),  VEL  (M) , DTVAL,  ErS,  DVSaPE,  HTl,AX 


1 
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1 


Si'b^PUT  Jfifc  St  Tt' A T (MATER,  TASS,  GAMMA,  RHOC"N,  N) 


3 

4 

5 

6 
7 
6 
9 

t 0 
1 1 
12 
1 3 
U 
15 
lb 
17 

1 A 
t 9 
20 
21 
?2 
23 

2 - 
25 

2b 

2t 

28 

29 

Jo 

31 

32 

33 

34 

35 
3b 

37 

38 

39 
«0 
“1 
«2 
43 

<JU 

45 

«b 

«r 

J8 

UQ 

50 

S| 

52 

53 
5a 

55 
5b 
57 

56 

59 

60 
fri 
62 


C 

C 

c 

c 

c 

c 

c 

c 

f 

c 

c 

c 

c 

c 


c 

c 

c 


10") 


1 o02 
1003 


1 o c a 


C 


mo 


c 

: 

c 

c 

c 

c 

c 

c 


this  routine  p»"vn.ES  a vectorized  double  precision  equation 

OF  STA7E  CALCULATION  rCR  a r I N C • WHEN  the  EQUATION  OP  state  is 

cmangfd,  a number  of  othep  routines  rust  re  edified  as  hell. 


EQUATION  OF 

HATFFfl) 


RASSC(I) 

GAhmaC(I) 

eftccd 

RHOCCT) 


STATE  ...  PH  0 s MHOC  ♦ CFPE/ENTCJ»*(1/GAHPAC). 

* cell  IDENTIFIER  = L.MM  HHEPE  0 < L < 10  IS  TH£ 

Layer  number  and  0 < NR  < 100  IS 

THE  MATERIAL  IDENTIFIER. 

* CELL  MASS  * PHfi (I)*LAH(I) 

* CELL  ADIABATIC  GAS  CONSTANT  - FIXED  DURING  ADINC 

= cell  ENTROPY  - CONSTANT  DURING  ADINC  HYDRODYNAMICS 
s DENSITY  CONSTANT  IN  THE  EQUATION  OF  ST aTE 


FARANETER  Npt  s 20? 

F E A l « 6 DCLR(NPT),  PELV(NPT), 
PEAl*P  maTLR(NPT),  MASS (NPT  ) , 

H E AL  *6  MASSC(NPT),  ENTC(NPT), 

REAL»8  Gamma  c (nPT  ) , 

COMMON  /ADICOM/  matERC.  MaSSC 


C Rm  A X , D P M J M 

GAMMA (NPT) , RhOCON(nPT) 
RHOC(NPT),  MaTEPC(NRT) 
DElRHO (NPT) 

GammaC  , E N T C , RHOC,  NCELLS 


CHECK  THE  INPUT  TO  SETmaT  FOR  REASONABLENESS. 
NPTO  S NRT 


IF 

( N . LC • 2 .OR.  N.GT. 

mPT-2) 

HPITE  (b, 

1001)  N 

, NPTO 

IF 

(N.lE.2  . Op . N.GT. 

nPT-2) 

STOP 

Call  maxmin  rATER(2) 

. N,  CRM  Ax , 

n ax, 

DRM I N , 

IhIN) 

IF 

(CRMIN  .LE.  R.OCO) 

HRITE 

(6, 

1002) 

DRW I N , 

IM  IN 

IF 

( r R M I n .LE.  P.OIP) 

stop 

CALL  maxmin  ( m ASS  ( 2 ) , 

N,  DRMAX, 

IMAX, 

r P M i N , 

I MIN) 

IF 

( D R M I N .LF  . L 0 ) 

write 

(6, 

1003) 

DPMJN, 

I M I N 

IF 

(DRI1IN  ,LE  . L.t'CO) 

STOP 

CALL  maxmin  ( R h 0 C 0 N ( 2 ) , N, 

DRMAX 

, I MA  X 

, drmin 

, ININ) 

IF 

(CRMIN  .LT.  O.L'DP) 

WRITE 

(6, 

10  0 4) 

DRM 1 N , 

I M I N 

IF 

(DRNIN  .LT.  O.OCO) 

STOP 

format  C(>SETuaT  input  PR051EM.  N out  of  Range.  I“,  14* 

' CALCULATION  STOPPED.') 

format  Pi>5ET*AT  INPUT  PROBLEM,  CELL  MATERIAL  NEGATIVE 
lPri2.4,  ' at  CELL  '*  14.  ' CALCULATION  STOPPED.') 

fcpmat  ccsetmat  input  problem,  cell  mass  negative 


1PD12.4,  ' AT  CELL  ',  I", 


CALCULATION  STOPPED. ' ) 


format  CoSETMaT  input  problem.  CELL  RHOCON  negative 


1PC12.4,  ' AT  CELL 


iu, 


calculation  STOPPED.  ' ) 


NP  s N * 1 
DO  ICO  I « 2,  NP 

natERC(I)  = NATER(I) 
MASSC(I)  * m A 5 s ( I ) 
GAMMAC ( I ) * CAMHACI) 
FHCC(I)  « RHCCON(I) 
NCELLS  ■ N 
RETURN 


ENTRY  SFTEOS  (RHP,  PRE,  N) 


SETEOS  COMPUTES  THE  ENTPOPY  CELL  CONSTANTS  (ENTC)  GIVEN  KNO*N 
VALUES  OF  THE  DENSITY  and  PRESSURE  IN  THE  CELLS.  SETEOS  IS  USED 
at  the  BEGINNING  Of  CALCULATIONS  FOR  INITIALIZATION  AND  DURING 
CALCULATIONS  WHENEVER  NON-IDEAL,  SOURCE,  OR  DISSIPATIVE  effects 
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63 

64 

65 

6 o 

67 

68 
69 
TO 

7 I 
’2 
73 
7 u 

75 

76 
7 7 

78 

79 

80 

91 

92 
83 

R« 

85 

86 
87 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


1005 


HAVE  CHANGED  The  CELL  ENTROPIES.  SINCE  aDINC  IS  PREDICATED  ON  THE 
CONSTANCY  of  The  LaGRANGIaN  CELL  ENTROPY  DURING  A FLUID  TIMESTEP. 
THIS  AMOUNTS  to  A form  of  TIMESTEP  splitting. 

R H 3 ( I ) S DENSITY  of  MATERIAL  IN  CELL  I (I  * 2,  N+1)  (GM/CC) 

THFSE  QUANTITIES  ARE  GIVEN  AS  INPUTS  TO  SETeOS 
PPEU1  * PRESSURE  GIVEN  IN  CELL  I BETWEEN  INTERFACES  I AND 

1-1.  THESE  VALUES  ARE  INFUT  TO  SETEOS  (ERG/CC) 

N s Num6ER  OF  INTERIOR  CELLS  IN  THE  SYSTEM 

RE  AL  *8  RHO(NPT),  PRE(NPT) 


CHECK  THE  INPUT  TO  SETEOS  FOR  REASONABLENESS. 

NPTO  = NPT 

IF  (N.LE.2  .OR.  N.GT.nPT-2)  WRITE  .(6,  1005)  N,  NPTfl 
IF  (N.LE.2  .OR.  N.GT.mPT-2)  STOP 

CALL  MAXMIN  (PH0(2),  N,  E'RMaX,  IMAX,  DPMIN,  IMIN) 

IF  CDRMIN  .LE.  O.OCO)  WRITE  (6,  1006)  ORPIN,  IMIN 
IF  ( D R M I N .LE.  0.000)  STOP 

CALL  MAXMIN  c PRE  c 2 ) # N,  DRMAX,  IMAX,  DRM I N , IMIN) 

IF  (ijRHlN  .LE.  O.CDO)  WRITE  (6,  1007)  ORMIN,  IMIN 
IF  ( D R M I N ;LE . O.OCO)  STOP 

FORMAT  C0SETE8S  INPUT  PROBLEM,  N OUT  OF  RANGE.  I«,  l»> 

I 1 CALCULATION  STOPPED. ' ) 

FORMAT  C'OSETEPS  INPUT  PPOELEM.  CFLL  DENSITY  NEGATIVE  '/ 

1 PD  1 2 . R i 1 AT  CELL  ',  I«*  1 CALCULATION  STOPPED.') 

fppmat  ('"Setfcs  input  problem,  cell  pressure  negative 


1PD12.4,  ' at  CELL 


I«» 


CALCULATION  STOPPED.  ' ) 


NP  ■ M ♦ 1 

DO  ?On  I * 2f  NP 

OELRHO(I)  = DN'AXI  (1.0D-JO,  (RHP(I)  - RmOC(I))) 
CALL  DUBLOG  tDELPH0(2),  ENTCC2).  N) 

CO  ? jo  I * 2,  NP 

ENTC  Cl)  = LNTC(I)*GAMMACCI) 

CO  215  I = 2,  NP 

ENTCCI)  = DEXP(ENTCd)) 

O'1  220  I = 2,  NP 
ENTC(T)  = prl  C I ) /Fl.TC  f I ) 

RFTUPN 


ENTRY  IJSEFOS  [RHO,  PPE,  LAMfOS,  DLAMUP,  N) 


USEEOS  TAKES  The  Set  of  CFLL  PRESSURES  (PRE)  and  CALCULATES  the 
CELL  DENSITIES  EXPECTEC  FOP  that  PPESSURE  based  on  The  EQUATION 
0 F STATE  CONSTANTS  FOP  THE  CELL  STORED  IN  COMMON  B|_CCK  /ADIC0M/. 
the  EXPECTED  CELL  VOLUME  (LAMEOS)  IS  COMPUTED  FOR  E*CH  CELL  AS  IS 
THE  JACOBEAN  DERIVATIVE  Dlamdp. 


PHOCI) 
PPE(I) 
LAuEOS(I) 
CL  a M C P C I ) 


= CENSITY  Of  MATERIAL  IN  CELL  I (I  = 2,  N*1)  (GM/cC) 
THES1  QUANTITIES  ARE  COMPUTED  AS  OUTPUTS  OF  USEEOS 
a PRESSURE  GIVEN  IN  CELL  I BFTWEEN  INTERFACES  I AMD 

i-t.  these  values  are  input  to  useeos  cfpg/cc) 

* Tnr  CELL  VOLUME  REQUIRED  BY  THE  EQUATION  OF  STATE 
GIVEN  the  cell  PRESSURE  PRE(I)  and  e.p.s.  constants 

* Rate  of  change  of  CELL  volume  with  PRESSURE  In  the 
EQUAtION  OF  state 

* MUM(,EP  TF  INTERIOR  CELLS  IN  THE  SYSTEM 


■BMW 
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r 


i 


125  RFAL*8  U*EnS(NPT),  DLAPPPEMPT) 

126  C 
' 27  C 

12B  C CHECK  THE  INPUT  T°  USE t T S FBP  PE ASBN ABLE NESS . 

129 
1 JO 
1 Jl 
132 
1 JJ 
t 3 « 

135  loop 

13fc  I 

137  1 u o 9 

1 3P  1 

139  C 
HO 
1 “1 

1«2  300 

' 9 3 
t«u 

195  320 

1"6 

197  390 

tup  C 


1 UO 

C 

evaluate  thf 

EQUATION  or  STATE  RPO  s ph«C  ♦ (P9E/F-NTC ) ** ' /GaH”AC 

1 5C 

D1  360  I 

* 2#  M 

1 5 1 

KHO  f I ) s 

Pt"*C  1 1 ) ♦ DELI  PCI) 

152 

360 

L *"E  r S C I ) 

* iAssctn/pHflcn 

153 

DO  3S0  I 

* 2,  M 

I5U 

380 

Dl A PDp C I ) 

s - LAHEBS(I)*CELr(I)/(rP8(I)*ppE(I)«GAM!1AC(I)) 

155 

RE'  TURN 

156 

END 

NPTP  s *-ip  T 

IP  (N.Lf.2  .OP.  r..GT.K'PT-2)  WRITE  ( 6 , 100B)  t.,  ‘IPTP 
ip  (f.lE.2  ,0«.  N.GT.mPT-2)  STOP 

CALI  W*  A X M 1 14  ( PRF  ( 2 ) , N / DPP  AX.  IHAX,  DP”!*,  I^IN) 

IP  (DPMI'1  .It  . 0.010)  WRITE  ( 6 , 1009)  pRt'Tl.,  JMIN 
IP  ( r p f*  I f I .Lt.  O.ODO)  STOP 

FORMAT  ('PlJStECS  I fiPUT  PPOPLCM,  t.  CIlT  OF  PANC-E.  I«,  19, 

» C^LCULATlOf.  STOPPED.  ' ) 

FORMAT  ('OUStEPS  INPUT  PROniEM.  CELL  PPFSSURE  t.E  C-  A T I VE  >, 

1 PE  1 2 . 9 , ' AT  CELL  '»  19#  ' CALCULATION  STEPPED.’) 

M u N ♦ 1 
DO  3C0  I * 2,  N1 

DELV(I)  * Dt,AXi(PPftI),0.000)/LNTr(I) 

CALI  DUPLOG  (0ELV(2),  PELP(2),  N) 

DO  320  I s 2#  Ut 

DFLP(I)  = rPLp(I)/GA"MACCn 

0»  3«n  I * 2#  M 

L P Lp  C I ) s CEXptOELECn) 


E 
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V 


1 c Al'INC  - SIMPLIFIED  FIVE  L A YEP  TF$T  PROGRAM  march  1979 

2 C 

3 C AFIAPATIC  THPOUGH  INCOMPRESSIBLE  RL"w  I»  O'.E  DIMENSION  Sa^S  SHOO 

9 C 

5 C This  PPSC-Pap  DE“ONSTPaTES  USE  BF  The  Ari^C  PACKAGE  TB  S8LVE 

6 c A NL'HBEP  BF  8NE"Dl  ''E  NS  1 0‘  AL  , L AGP  ANGI AE  FLUID  DYNAMICS  pP®BLEMS  . 

7 C A i-ir.r  CLASS  BF  PRCBLE^S  can  fl£  SET  up  USING  THE  data  PPIV£N  IN- 

0 C ITIALJZFP.  ELEXI6LC  DIAGNOSTICS  EBP  EACH  Bp  T H£  FLUID  LAYERS  »FE 

9 c included.  special  attention  has  e e e n paip  tb  treating  dfnsjty  dis 

10  C CONTINUITIES  AND  ZONE  SIZE  r I sc BNT I NU I T IF S ACCURATELY,  the  PPESEN 

11  C VFrSION  nr  ADI'  C and  its  TEST  PR"GPA"  ARE  WRITTEN  E*!T7RELy  In  faU 

1?  C BIT  (DOUBLE  PRECISION)  FLOATING  pbinT  ARITHMETIC,  THUS  CONVERGED 

1 3 c t p bettff  than  i part  if  io**t  is  possible  for  problems  *Ith  near 

14  C INCBhRPESSIdIL1TY  aND/CP  EXTREME  DENSITY  DISCONTINUITIES  + hIcm 

15  C PEQUJ RE  THIS  ACCllPACY.  for  SCMEWm4t  LESS  EXTREHE  FLUID  SYSTEMS, 

16  C 32  RIT  RLBATING  PC I NT  COMPUTATIONS  S lOULO  RE  ADFQi'aTE  AND  T-f 

17  C C "RPE SPONDI NG  SINGLE  PRECISION  TRIPIaGBNAL  S"LVEPS  ARE  AVAILABLE. 

ie  c 

19  C ADINC  has  BEEN  CONSTRUCTED  AS  a UTILITY  RaCKAGE  TO  ADVANCE  TH 

20  C FOUR  HYPFBnVllAMIC  VARIAPLRS.. 

?1  c. 

22  C R»D(I)  = POSITION  (RADIUS)  BF  The  I-TH  CELL  INTERFACE  (Cm) 

23  C VELdl  s VELBCITY  or  thf  i_tm  CELL  INTERFACE  (C^/SeC) 

24  C RHO(l)  * DENSITY  IN  CELL  I BETWEEN  INTERFACES  1,1+1  (Gm/CC) 

25  C FFe(I)  * PRESSURE  In  Thf  i-TH  c BHRUTa t I Oral  CELL  (tRG/CC) 

?o  C 

2-  C LAGpANGlAN  fluid  DYNAMICS  EQUATIONS  ARE  SOuvED  INCLUDING  a RlEX- 

28  C I BlF  E 0 U A T I B t , OF  STATE  wt  ich  CAN  VARY  FROM  CELL  TO  CELL  IN  The 

2R  C ClSCRETIZEr  REPRESENTATION  BF  THE  FLUID.  ADP'C  has  BEEN  C‘ST  I N T f 

30  C A FORM  RESEMBLING  That  wf  An  opdINary  DIFFERENTIAL  ECUaTiBnS 

3 1 C PACKAGE.  ThE  USER  CAN  REQUEST  INTEGRATE  T"  a CFRTAIn  Tjme  jnde- 

3e  c ffmoant  bf  the  number  bf  ctclfs  rfduired  in  adinc.  the  usfr  alsb 

33  C HAS  CONTROL  BF  VAUIBUS  ERROR  AND  INTEGRATION  RaRawETE=S  WITHOUT 

34  C HAVING  TO  PLUNGE  INTO  ThE  pOWFLS  BF  THE  SBeutIBN  m£TM"D  ItSELF, 

35  C THF.  EQUATIONS  SBlVED  »r.f 

3o  C 

37  C ii(RAD)  D(VEL)  - ! 

30  C = VfcL,  s GRAD  (PrE), 

39  C OT  CT  RhB 

4Q  C 

4(  C A No  TH£  EGUAT I BN  wr  state  .. 

«2  C 

43  c ( PRF  ) **  1/GAMHAc 

44  C RHJ  x PHOC  ♦ ( ) 

95  C ( Eetc). 

«6  C 

97  c THE  T(  ST  °FBgRAn  is  ARRANGED  TB  HANDLE  UP  TO  RJVE  DISTINCT 

90  c LAYERS  bf  FLUID  cBMPOSEC  OF  up  TB  ?C0  INDIVIDUAL  FINITE  dIFPErEHC1 

hr  c celis.  facn  cell  is  lagrangian  and  hAs  several  quantities  that  Ari 

5.)  C CBnSePVFD  MOVING  +ITH  the  FLUID  AS  LONG  AS  DIFFUSIVE  AND  BthEB 


75 


51 

52 

53 
5u 
55 
5b 

57 

58 

59 

60 
61 
62 
63 

6i| 

65 

66 

67 

68 

69 

70 

71 

72 

73 
79 
75 
7 1 

77 

78 


83 

6a 

85 

86 

87 

88 
89 
«0 

<»2 

03 

9a 

05 

96 

97 

98 

99 
100 
101 
102 
103 
toa 

105 

106 

107 

108 
106 
HO 
1 1 : 
1'2 


NOn-IDEAL  and  EXTERNAL  SOURCE  EFFECTS  ARE  IGNORED.  THESE  CONSTANTS 
ARE  COMMUNICATED  THROUGHOUT  THE  ADINC  PACKAGE  IN  common  BtOCK 
/ADICOM/.  THESE  "CONSTANTS"  VARY  FROM  CELL  TO  CELL  ACCORDING  TO 

the  initial  conditions,  following  are  the  definitions  of  these 

QUANTITIES.. 


materc(I) 


masscci) 

GAHMAC(I) 

ENTCCI) 

RHOC(I) 

NCELLS 


■ CELL  IDENTIFIER  * L.MM  WHERE  0 < L < 10  IS  THE 

LAVER  number  and  0 < mm  < 100  IS 
the  MATEPIAL  IDENTIFIER. 

“ CELL  MASS  ■ RH"(I)*LAM(I)  - HELP  CONSTANT  IN  aDINC 

* cell  aciabatic  GAS  constant  - HELP  fixed  in  ADINC 

* cell  ENTROPY  - CONSTANT  DURING  aDINC  HYDRODYNAMICS 
« DENSITY  CONSTANT  IN  THE  EQUATION  fF  STATE 

e NUMBER  of  CELLS  "F  FLUID  IN  THE  CALCULATION 


VARIABLE  gamma  and  ENTROPY  aRE  USED  IN  EACH  LAGRANGIaN  CELL 
ANC  AN  IMPLICIT  PRESSURE  ITERATION  ensures  LINEAR  STABILITY,  IP 
NOT  ACCURACY,  FOR  TIMESTEPS  LONGER  than  the  COURAN'T  CONDITION,  the 
NOnLINEaP  TERMS  A«E  ITERATED  with  a QUADRATICALLY  CONVERGENT  ALGO- 
RITHM, ADINC  is  DESCRIBED  IN  DCL  MEMORANDUM  REPORT  « , 

1979.  the  TPIDIAGONAL  SOLVERS  USED  ATE  DOCUMENTED  IN  NPL  memoran- 
dum REPORT  “3908,  NOVEMBER  1976. 

PROBLEMS  IN  ONE  OF  FOUR  GEOMETPIeS  CAN  BE  SFT  up  FOR  “DlNC  BY 
CHANGING  THE  INTEGER  ALPHA  IN  THE  CALL  TO  SETGEO.  ALL  NON-IDEAL 
EFFECTS  HAVE  TO  BE  TREATED  SEPARATELY  EITHER  BY  IMBEDDING  8R  BY 

timestep  splitting. 


79 

C 

alpha 

s 1 

CARTESIAN  COORDINATES. 

«C 

C 

alpha 

= 2 

cylindrical  COORDINATES. 

81 

C 

ALPHA 

* 3 

SPHERICAL  COORDINATES. 

F2 

C 

alpha 

* a 

POWER  SERIES  (NOZZLE)  COORDINATES 

the  boundary  Conditions  treated  ry  apikc  ape  quite  general. 
the  POSITIONS  and  VELOCITIES  OF  the  region  BOUNDING  INTERFACES 
RAD ( 1 ) AMD  R A D ( N 1 1 ) CAN  Rf  EXTEPNALLY  DETERMINED  FUNCTIONS  OF  time 
and  OTHER  PHYSICAL  VARIABLES  DURING  THE  CALCULATION.  THE  BOUNDARY 

conditions  are  communicated  to  the  adinc  package  via  thf  f?up 
arguments. . . 


RRnE  w 
RLnE*. 
VRnE  W 
VLNEw 


* RIGHT  BOUNDARY  POSITION  RADCN+l)  AT  END  Of  TI“ESTEP 
= LEFT  BOUNDARY  POSITION  p a D ( 1 5 AT  OF  TIMtSTFP 

s RIGHT  BOUNDARY  VELOCITY  vFLCM+t)  AT  END  Of  T 1 1 ; E b T E P 
“ LEFT  BOUNDARY  VELOCITY  VEL(l)  at  END  Of  timestep 


SEVERAL  AUXILIARY  VARIABLES  ARE  USED  °Y  ADINC  and  The  CONTROL- 
LING TEST  PROGRAM  WHICH  ShOLLD  ALSO  BE  EXPLAINED  TO  The  user. 


N s NCELLS 
DElTaT 


ISTEP 

EPSRO 

EpSvO 

MDaMP 

EPS« 

EPSV 

GEOMCOd  ..5) 

APEA(I) 

PADC(I) 


NUMDEP  3E  FLUID  CELLS  IN  The  ADINC  INTfGPatION 
TiMt  INTERVAL  for  the  adinc  INTEGRATION  Whlcw  may 
SUBCYCLE  UP  TO  100  times  INTERNALLY  IE  NEEDED  FOR 
ACCURACY  op  STABILITY. 

TIMESTEP  NUMBER  USED  by  adinc  FOR  IDENTIFICATION 
EXpLICITNESS  PaRaMFTER  FOR  THE  POSITION  INTEGRATION 
EXPLICITNESS  PARAMETER  FOR  the  VELOCITY  INTEGRATION 
NuMRFP  OF  CYCLES  AT  THE  BEGINNING  OF  a CALCULATION 
IN  WHICH  ADDITIONAL  DAMPING/SMOOTHING  IS  APPLIED. 

fxplicitness  parameter  fop  rad  last  used  bv  adinc 
explicitness  parameter  for  vel  last  used  bt  adinc 

array  OF  FIVE  COORDINATE  COEFFICIENTS  for  ALPHA  c U 
area  of  the  1-th  CEIL  INTERFaCF 
position  of  the  i-th  cell  centep 
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U5 
1 1 a 
115 
lit 
117 


131 

132 

133 
139 
135 
1 3t 

137 

138 

139 
tuo 
i«i 
1«2 
1«3 
1UU 


LAMCCI)  S VOLUME  0 F THE  I-TH  finite  DIFFFBE'CC  cell 

EKIN'ECI)  * KINETIC  energy  DENSITY  in  THF  I-tn  cell 

EThRM(I)  ■ THERMAL  ENERGY  DENSITY  IN  Thf  1-th  CELL 


1 18 

PARAMETER 

NPT  S 202 

1 1 9 

INTEGER 

alpha 

120 

LOGICAL 

LTPPT, 

LCHEM, 

ltcnd. 

L CIFF 

121 

LOGICAL 

LZOnE 

122 

RFALoR 

LANPCMPT), 

DABSV(NPT), 

ETOTL(NPT) , 

THFRACCNPT) 

123 

REAL*8 

FTHPM(NPT) , 

FKINF(NPT), 

EPSP, 

F P S V 

I 2U 

RE  AL  *8 

dtmax. 

DTMIN, 

PELTAT, 

time 

125 

REAL*8 

P L N E w , 

VLNEW, 

RRNEw, 

VRk'E  h 

12b 

RFAL*8 

EPSRO, 

FPSVO, 

DTIME, 

T I ME  n ( 1 1) 

127 

REAL  *6 

PARC* 

dtprf, 

dtval. 

G E 0 M C 0 C 5 ) 

126 

REAL  *6 

PAD(NPT)  , 

VEL  (NFTJ  , 

RHO (NFT ) , 

P R F C N P T ) 

129 

REAL*8 

AREA  CNPT)  , 

PADCCNPT) , 

LAMC (NPT ) 

130 

real  m 

POLDtNPT), 

VOLDCNPT), 

PAVGCNPT), 

vavgcnpt) 

the  FOLLOWING  DECLARATIONS  APPEAR  IN  ADINC  AND  THF  E.o.S.  ROUTINES 
FOR  THE  COMMON  BLOCK  ADICOM  WHICH  CONVEYS  AND  STORES  E.O.S.  INFOR- 
MATION. it  is  needed  in  the  main  program  for  diagnostic  reasons. 

realms  nASSC(NPT),  ENTC(NPT),  RHOCCliPT),  haTERCCnPT) 

REAL*8  GaMHAC(NPT) 

COMMON  /AOICOM/  MaTERC.  MASSC,  GaMMac,  E n T C , RhOC,  ncells 

DATA  MAxSTp,  IPRIAJT,  IPLOT,  ALPHA  /26,  1,  0,  1/ 

DATA  PTM I N , D Tn  A X , N /l.OD-1,  1.00-1,  10/ 

DATA  EPSRO,  EPSVO  /0.5CDO,  0.50DO/,  TImEC  /HaO.OOO/ 

DATA  GEONCO  /t.OfO,  u.OOO,  P.OCO,  O.ODO,  O.ODCV 

data  lzone  /.falsf./,  lchem  /.false./,  ldiff  /.false./ 

data  ltcnd  /.false./,  ltppt  /.true./,  noa<jp  n / 


195 

19b 

l«7 

i«e 

199 

150 

151 

152 

153 
159 

155 
15t 
157 

156 
159 
) 60 
161 
162 
163 

1 69 

165 

166 

167 

168 
1 69 

170 

171 
>72 
173 

ITU 


NAMELIST  /CANTRL/  MAXSTP,  IPRINT,  ALPHA,  DTMIN,  DTMAX,  n, 

1 EPSPO,  EPSVO,  GEOACC,  I.ZONF,  LCHFm, 

2 ItlFF,  LTCND,  LTPRT,  mdamP 
C 

1000  FORMAT  CM  aETER  step  NO.  '*15,'  DT  s Dlfc.P, 

1 ' THE  TIME  IS  ',  C 1 6 . P , /) 

1001  FORMAT  C 2 X , I 3 , 1P9D12.9,  0P2F8.9) 

1002  FORMAT  C I POSITION  VELOCITY  IENSITY  i, 

1 ' pressure  ENERGY  AVG  PRESS  AVG  VELOC  '* 

? ' CELL  VOLUME  fNTROFY  THFPaC  Gamma  /) 

1003  FORMAT  C'CNAMELIST  /CONTRL/  DEFAULTS  ...'} 

1 0 0 9 FORMAT  C'ONAMFLIST  /CONTRL/  UPDATES  ...') 

1005  FORMAT  ('  at  STEP  '.  15,  ' THE  RFOUIPfO  TIhfStEP  1 P C 1 2 . 9 , 

1 ' BECAME  smaller  than  the  minimum  step  ',  P l 2 • 9 ) 

?00 1 FORMAT  COTIMlNGS  at  CYCLE  15,  ' TI“E  s iDDlS.9,  /, 

1 ' RfcZONF  TIMESTEP  CEIL  PRINT  l avEc  PP^'T 

? i GRAPHICS  CHFMISTRY  DIFFUSION  CONDUCTION 

3 ' ADINC  ') 

2002  FORMAT  CSX,  1P10D12.9) 

?003  FORmAT  C TOO  many  REZONES  at  cycle  19,  ' NX  s ' , Iu) 

c 

c 

c ********************************* 

C CPnTROL  PARAMETERS  ARF.  INITIALIZE^,  change  F"p  DIFFERENT  CaSFS. 

c a*******************************  * 

call  common 

TIME  ■ O.eto 
WRITE  (6,  1 U C 3 ) 

HR  I Tc  (6,  CPN'TfcL) 
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175 
17fe 
177 
1 78 
1 7q 
1 80 
181 
tea 
183 

18U 

165 

ISfc 

187 

188 
1 «<? 
1 90 
I’l 
192 
19J 

19a 
1 *»5 

196 

197 

198 
l 9 q 
200 
201 
20? 
?03 
20a 
205 
20fe 

207 

208 
2no 
210 
211 
212 
213 
210 

215 

216 

217 

218 
216 
22  0 
221 
222 
223 
229 

225 

226 
2?7 
226 

229 

230 

231 

232 

233 
23a 
235 


FEAN  (5*  COt-TPl.) 

WRITE  (b#  190“) 

WRITE  (6,  C 0 N T P l ) 

L'TPPE  = OSQRT  (OTMIN*r>TMAX) 


initialize  the  physical  data  fop  the  calculation  including  the 
GRIP#  MATERIAL,  a 'iri  FLUID  PROPERTIES.  a USep  may  wish  t«  DEVELC'Pc 
His  own  rpfFLEH  IriITIALIZFR. 

**************  * ****************** 

call  I r.  I T A L CmaTEFC,  GAM*'ac,  PHOC , RAD,  VFL  # phO,  PPE#  nj 
UP  S N * 1 
RENEW  » R*0(1) 

VLNFW  s VEL(l) 

KRNFW  s RAL(llP) 

V P N F W s VCL(NF) 

A*********************#,*****,**., 

INITIALIZE  the  ADlNC  GE8MFTPY  AND  F QUA T 1 ON.WF-ST ATE  MODULES  AND 
THE  VARIARLE  EXPLICITNESS  PAPAHETEKS  within  ADINC.  whfn  HATEPC, 
MA$SC,  GAMMAC,  AND  RHOC  ARE  DEFINED  IN  PLACE  AS  They  aPE  heR£,  Thf 
CALL  TP  SETNaT  BtLflW  IS,  STRICTLY  SPEAKING,  S"PF PFU'&US. 

A***************************  * * * « » 

CALL  SETGEO  (Alpha,  GFPhcO) 

CALL  USEGEO  (Rad,  area,  FADE,  LAMC,  N) 

DP  10  I s 2 , Np 
VPlDCII  b V e L C I j 

10  MASSC(I)  e RHfl(I)*LAHC(IJ 
VPLDCU  B vel  c 1 3 

CALL  SETHAT  (HaTERC,  HASSC,  GAMMAC,  RHPC,  N) 
call  SEtEDS  (RhP,  PRE,  N) 

CALL  SC 7 EpS  (EPSRO,  EPSVO,  MCAMP,  EP5R,  EPSV) 

hi******************************* 

THE  LPPP  PVER  TIMESTEPS.  THE  VaRIPUS  SUBSECTIONS  ARE  TIMED  FpR 
DIaGHOSTIC  PURFOSES. 

********************************* 
DP  9999  ISTEP  s t,  MAXSTP 

WARNING  ...  THE  FOLLOWING  BLOCK  PF  CODE  C8NTRPLS  THE  PRINTPUT 
INTERVAL  AND  the  EXPLICITNESS  PARAMETERS  near  TURNaRPUND  In  the 
LINUS  SIMULATION.  THEY  are  PROBLEM  SPECIFIC  AND  SHPULD  BE  PEHPVED 
OR  BYPASSED  FOR  GENERAL  APPLICATIONS. 

IF  (ISTEP. GE. 965  .AND’.  I STEP . L E . 7 1 0 ) 
l CALL • $ETEPS  (O.ODO,  O.ODO,  MdAMP,  £PSP,  EPSV) 

IP  (ISTEP. EQ. 950)  IPRINT  ■ 5 
IF  (ISTEP. EQ. 9b5)  IPRJNT  s 1 
IP  (ISTEP. EO. 975)  IPRINT  b 5 
IF  (ISTEP. EQ. 7icj  IPRINT  s 50 

********************************* 

LAGRANGIAN  REZONE  FACILITY  »«  TO  BE  DEVELOPED 
********************************* 
IP  (.not.  lZONE)  GO  TO  15 
CALL  SECOND  (1,  DTIME) 

* * * A 
AAA* 

CALL  SECOND  (0,  DTIME) 

TIMED ( 1)  = TlMED(  1)  ♦ DTIME 
15  CONTINUE 
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V 


236 

237 
23R 
239 
2«0 
2«1 
2«2 
2«J 
2«4 
2«5 
2<*6 

247 

248 

249 

250 

25 1 

252 

253 

254 

255 
25b 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 

267 

268 

269 

270 
2?1 

272 

273 

274 

275 

276 

277 

278 

279 

280 
2"1 
262 
263 

284 

285 

2«b 

287 

208 

209 

290 

291 

292 

293 

294 

295 

296 

297 


C ********************************* 

C CALCULATE  THE  TIMESTEP  BASED  ON  FLOW  VELOCITY  AND  EXTERNAL  LIMITS. 

C A**************************.**.** 

CALL  SECOND  (I,  DTIME) 

DP  20  I ■ 1»  Np 
DABSV(I)  ■ 0a8S(VEL(I)J 

20  DABSV(I)  ■ DNAX1  CDABSV (13,  DABS ( VOLD ( I ) ) ) 

CALL  DTFlPW  (RaD,  UABSV,  DTVAL,  N) 

DELTAT  » DMlNt  CDTMAX,  0.35D0*DTVAL) 

DELtAT  ■ PM IN ' (DELTAT,  1.1D0*DTPRE) 

IF  (DELTAT  ,lT.  DTt'IN)  WRITE  (6,  1005)  ISTEP,  DELTAT,  OTMIN 
DEL7AT  * DMAX1  (DELTAT,  DTMIN) 

DTPRE  a OELTAT 

CALL  SECOND  (<’,  DTIME) 

T 1 HEP ( 2)  s TIMED  t 2)  + DTIME 
C 

C <L«4  ******************************* 

C , THE  RESULTS  and  DIAGNOSTICS  are  PRINTED  by  CELL  and  by  LAVE».  the 
C UTILITY  EPGPRT  USES  the  CELL  IDENTIFIER  MATERC(I)  to  distinguish 
c layer  BOUNDARIES  Fqp  certain  IMPORTANT  CONSERVATION  S"MS. 

c ********************************* 

IF  CIPRINT  .fU,  0)  GO  TO  30 
IF  (MODdSTEP-l.IPRINT)  .NE.  C)  GO  TO  30 
C 

C *********8************ 

C CALCULATE  THE  ENERGIES  AND  OTHER  QUANTITIES 

C ********************** 

CALL  SECOND  (1,  DTIME) 

CALL  USEGEO  (RaO,  ARCa,  RADC,  LAMC,  N) 

DO  25  I ■ 2,  Np 

ETHRM(I)  * PEE  (I)  / ( G ANN  AC  ( I ) - l.t'DO) 

C 

C THE  FOLLOWING  FUNCTION  EVALUATES  THE  ENERGY  INTEGRAL  POP  A Ga*m* 

C EQUALS  1/2  "SOLID"  EQUATION  OF  STATE. 

IF  (RHPC(I)  ,NE.  O.ODO) 

1 C»PG  * r-HAxi  (RHO  ( I ) /RHPC  ( I ) , 1 .OOoOCOf'OOOOlD't) 

IF  (GANNAC(I)  ,EQ.  0.5DO)  ETHRM(I)  s -PRt  ( I ) * ( 1 . Or  !> 

1 - OaRGaPaTaN  (DSGRT  (PaRG  - t .000)  ) /DSOPT  (DARG  - l.oPOT 

C 

LAMP  ( I ) a (RAD(I)  - RACC(I))*(LA,'C(I)  ♦ 0 .5D0*  (PArC  ( I ) * 

1 PAP(I-1))*(AREA(I)  - AREA (1*1 ) ) )/(PAP(I)  - RAO(I-I)) 

E k I N E ( I ) B 0.5oo*RHO (I  )*(VEL (I )**2*LAHP (I ) ♦ VEL(I-1)**2* 

1 (LAMC(I)  * LAmP(I)))/LAMC(I) 

ETOTL(I)  C EMNF(I)  * FTHRM(I) 

25  THFRAC  ( I ) a E'ThRM  ( I ) /F  TPTL  ( I ) 

JSTEP  a ISTEP  - ' 

WP I TE  (6,  l«'0O)  JSTEP,  DELTAT,  TIME 
WRITE  (6,  1002) 

WRITE  (6,  1 r»0 1 ) (I,  RADd),  VEL(I),  RHp(I),  PPE(I),  ETCTUI). 

1 PAVG(I),  VaVG(I),  LAMC(I),  ENTC(I),  THFRAC(I),  Gam"AC ( I ) , 

? 1=1/  NP) 

WRITE  (6,  2°C")  ISTEP,  Tt*E 
WRITE  (8,  ZOOS)  (TIMER  (JJ),  JJ  a 1,  9) 

C 

C ************************** 

C CHECK  ON  the  INTEGRATOR  PrPFPPMANCE  VIA  APINCO.  T ml 

C PREVENTS  THE  RESET  OF  THE  NC ALL  POINTER  IN  ADTNC  TO 

C SM0"THING  IS  Oil!  Y ENABLED  ON  THF  INITIAL  "aTA. 

C ************************** 

CALL  ADlllCO  (!,  ISTEP) 

WRITE  (6,  ICO') 

CALL  SECOND  (0,  PTIME) 


******* 

A°GUMfcN t t 

ZPPR.  THUS 

******* 


*********** 
AS  DIAGNOSTICS. 

*********** 
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?99 

300 
3o  l 
30* 
303 
30a 
305 
30b 
jo? 
30s 

30q 

3 ' 0 
3 t 1 
31* 
3 ’ 3 
3 1 4 
315 
3la 
317 
3 1 B 
319 
3*0 
3*1 
3** 
3*3 
32a 
3*5 
3*a 
3*7 
32a 
3*9 

3 3 0 
3 3 1 
33* 
333 
33a 

335 

336 

337 
33a 
339 
3«3 
3«1 
3«2 
3«3 
iaa 
305 
3«b 
3«7 
3«8 
3 J9 
350 
«i 
35* 
353 
35a 

355 

356 

357 
35? 

359 

360 


T I MEC  ( 3)  = TlMFOC  3)  + t TI«E 

c 

C *«***«»***«**********..*»....,,«», 

C EnePGV  anC  PThtR  1'IAGNPSTICS  SHFLL  "Y  SHELL.  THE  MATERIAL  CAN  VARY 

C HlTHp;  A SHELL  RUT  THE  SHfLLS  RUST  BE  ENCOUNTERED  IN  NURErICaL 

C SEQUENCE  FRO"  l T*  AT  HOST  5. 

C «.*»*»**««*«*«*****»*******»**»** 

CALL  SECOND  C * « LTIHE1 

call  fpgppt  (Rad,  vel,  area,  larc,  rho,  pre,  n,  ekine,  ethrm) 
call  SCCONt  (!',  DTIME5 
T I HE  0 ( a)  = T 1 1 ‘ E C1  ( a)  ♦ CTIHE 

c 

30  CONTINUE 
C 

c *»»**»*«*******»«* 

c PERFORM  GRAPHICS  **  TO  BE  DEVELOPED 

C ****************** 

ir  (IPLOT  .EL.  n)  GO  TO  35 
IF  (HSDdSTFr-i  , PLOT)  .HE.  0) 

LALL  SECOND  (1 , CTIHE) 
c * » * * 

c * * « * 

call  second  (0,  OTIRE) 

TIMED!  5)  s TIMED!  5)  ♦ DTIME 
35  CONTINUE 

C 

c ********************************* 

C INTEGRATE  the  CHEMICAL  kinetic  Rate  EQUATIONS  **  TP  be  DEVELOPED 

C ********************************* 

if  (.NPT.LCHEH  ) GO  T 9 Uc 
CAUL  SECOND  (1,  DT I me  5 

c * * * * 

c * * * * 

call  SECOND  P,  D T I me  ) 

TIMED!  t>)  = TIMED!  6)  ♦ CTIHE 
«c  continue 

c 

c ********************************* 

c PERFBPM  DIFFUSION  CALCULATIONS  **  TO  BE  DEVELOPED 

c ********************************* 

IF  C.NOT.LLIFF)  GO  T«  50 
CALL  SECOND  C,  DTIREi 

C * * * * 

c * * * * 

CALL  SEC ONC  CO,  DTIME) 

TIMED!  7)  = TIMED!  7)  ♦ L T I M E 
5c  CONTINUE 

C 

C ********************************* 

C.  PERFORM  CONDUCTION  AND  ENERGY  ADDITIONS  **  TO  BE  DEVELOPED 

C ********************************* 

IP  C.NOT.LTCND)  GO  TO  60 
CALL  SECOND  (1,  DTIME) 
c * * * * 

c * * * * 

call  SECOND  CO,  DTIME) 

TIMED!  6)  = TIMED!  6)  ♦ DTIME 
60  CONTINUE 

C 

C ********************************* 

C CALCULATE  THE  AUGMENTED  PRESSURE  AND  UPDATE  THE  ENTROPY  IF  NflN- 

c ideal  physics  is  Included’,  save  the  old  pressure  and  velocity. 

c *********************  ************ 


*************** 
*************** 
GO  TO  35 
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361 

362 

363 

364 

365 

366 

367 

368 

369 

370 
37t 

372 

373 
37u 
375 
37e 

377 

378 
37q 

38  0 

381 

382 
363 
38a 
3f>S 

38e 

387 

3*e 

389 

39(J 

391 

3^2 

393 

39a 

395 

3^6 

3*7 

398 

3°9 

aoo 


C 

C 

aoc 


1 no 


C 

c 

c 

c 

c 


70 


c 

c 

c 

c 

c 


1 1 0 


c 


9999 


D8  aoo  J s ?,  NP 
* * » » 

* » * * 

continue 

ip  (LCHEM  .OR.  LDIFF  ,8R.  LTC^D)  CALL  SETEBS  (RH8,  PRE,  N) 

C*  100  I s 2,  HP 
PflLD(I)  = PPL(I) 

VBLt>  C I ) = VEL  ( 1 3 

vbldcij  = veu  til 

********************************* 

integrate  a hydrodynamic  tikestep  bf  length  deltat  after  resetting 
the  BOUNDARY  CONDITIONS  T 8 the  FNO  BF  THE  TIMESTEP. 

A******************************** 

IF  C.noT.lTPPT)  GB  TB  76 
CALL  SECOND  (’,  CTIHE1 

V P N E R = VF  L ( Mp  3 

R PNE  8 = PANEL  + DELTATaVKNEW 

vlnew  = vel(i) 

RLNFK  = RLliEF  + DEL  T A T*  VLNE  h 

CALL  ADIMC  (K’Ao,  VEL,  RHB,  PRE,  N,  DELTAT,  ISTEP, 

RRNe«,  RLNEW,  VRNEL,  VLNEW) 

C-ALl.  SEC  81  D (<>,  C T I HE  3 
TI*’ED(  R 3 = TIMED!  93  ♦ DTIH£ 

CONTINUE 

* ******************************** 

CALCULATE  T«t  AVERAGE  VELOCITY  AND  PRESSURE  ACTUALLY  USED  EY  ADINC 
TB  ADVANCE  The  PBSJTIBNS  f RAD 3 and  VELBcITiES  ( VEL  3 For  DIAGNOSIS. 
********************************* 
CALL  SETePS  (E.FSRO,  EPSVO,  “DAMP,  £PSR,  EPSV) 

DO  no  I s ?,  nr 

V a VG { 1 3 = EpSR*VOl 1(13  ♦ (I.ODO  - EpSR  3 *VEL  ( 1 3 

R A VG ( 1 3 = EpS VaPBLn ( 1 3 ♦ (1.000  - EPSV)*PRE(I3 

V A VG ( 1 3 = EpSP*VBLD  ( 1 ) + (I.ODO  - EpSp)*VEL(13 

TIME  * TIME  + DELTAT 

continue 

STOP 

LNr. 
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3 

4 

5 

t> 

7 

8 
9 

'll 
1 1 
' 2 
'3 
i u 
15 
I b 
17 
1 8 

1 9 
20 
21 
?? 
?3 

?4 

25 

2 o 

27 

28 

29 

30 
Si 

32 

33 

34 

35 

50 
37 

36 

39 
<10 
“l 
u2 
43 
«4 
«5 
4 6 
47 

4 8 

40 

5n 

5 1 

52 

53 

54 

55 

5 6 

57 

58 

59 

60 
61 
62 


C 

C 

c 

c 

c 

r 

r 

c 

c 


1 0 o 6 
1 ,:07 
1 o 0 8. 

1 it  b 9 

1 

2(,06 

2 0 0 7 
2 oo6 
2 DO  9 

1 

2 

4(1  Ofe 
4007 

4rop 

4005 

4 0 0 9 

4o  t 0 
“Oil 
4012 
C 
C 
C 


subroutine  tr igprt  (»*M  , vNEw,  area,  lam,  putw,  pne«, 

r;x,  kK>,  EK*) 

t p cpf T takf.s  tnf  cell  amt  interface  data  ait  uses  the  material 
ir-fNTirifR  A p k A V ■' A T L A C TO  RETARULATF.  VaPIBUS  PHYSICAL  QUANTITIES 
t.Y  lAYEP  PaThEP  T"An  r Y CELL,  the  DIAGNOSTIC  IS  PARTICULARLY  use- 
mil  f*r  HfTEnjGEiirous  calculations,  only  five  layers  are  permitted 

C UF CE n TL  Y f;UT  THIS  COULD  EL  INCREASED  EASILY  BY  CHANGING  SOf'E  "E 
the  loop  LIMITS  aP'L  FORMATS.  THANKS  TO  eLLIOC  TENT  EOR  DEVELOPING 
THIS  ROUTINE  ESR  * N I N C . 

EAPANFTEP  NPT  S ?02 

tfalap  R'fkcnpt),  Piifw(r.PT),  knfw(npt),  enew(NPT) 

REALaP  Pam  (NPT),  VNfA(NPT),  APE  A (NPT) , LAM(NPT) 

REAL*®  XK  W{  ■»  c 3 1 3 * XINEW(U),  XENEW(li),  XGAH(ll) 

k F A L * 6 XMASS(ll),  XPHS(ll),  X P RE  ( 1 1 ) , XVOL(H) 

pFaL*P  XivAtdn,  XVEL(ll),  XAPEA(il),  xentoI) 

RE  Al  *8  TgAH,  tent,  TRhOC,  XRHOCdl) 

PFAL*B  TNASS,  tinew,  tknew,  TeN£W 

Rf  Al *6  Trho,  TPPE,  TVBL 

K F a L * 8 FASSC(NPT),  FNTC(NPT),  PHOC(t.PT),  MaTERC(NPT) 

r-  E A L * 8 GAUNAc  (NPT  ) 

COUPON  /ADICBM/  uATERC,  HASSC,  Gammac,  ENTc,  RHOC,  1 cells 

FORMAT ( 1 T M F P M a L EPG  ',  1P6D17.A) 

FORMAT ( 1 R I N( T i c ERG  ',  lPbDl7.fl) 
r 0 R M a T C ' TOTAL  EFG  ',  !PfcDl7.e) 

FORNATC..',  l3x,  ' TOTALS  SHELL  1 ', 

• shell  i shell  3 shell  4 shell  5 •) 

FORMATC'  INTERFACE  POSITION',  1P6D17.8) 

FORMAT ( ' INTERFACE  VELOCITY',  1P6D17.6) 

FORMATC'  INTERFACE  APfA  ’,  1P6P17.B) 

format c 3x,  'Layer  bouncary',  7 x,  'inter  i',  iox,  'inter  2>, 

1 0 X , 'INTER  3',  1 OX , 'INTER  4',  IOX,  ' INTEP  5',  lOX, 

' H ter  b ' ) 

FORMATC'  VOL  Avg  R|iO  ',  lPbD17.8) 

FORMATC'  VOL  Avg  PRE  ',  1 pfcD 1 7 , 8 ) 

FORMATC  LAYER  VOLUME’,  1P6D17.8) 

FORMATC'  LAYFR  HASS  '»  1P6C17.8) 

F 0 R M A T C ' VOL  4VG  GAM  lPbD17.8) 

FOrhatc'  VOL  *VC-  ENT  ',  1P6D17.8) 

FORMATC'  VOL  *VG  RHOC',  1P6D17.B) 

FORMAT  ( 5X , /) 


THE  VARIOUS  SUMMANCS  ARE  INITIALIZED  TO  ZERO. 
NP  s MX  ♦ 1 
T''ASS  = O.ODO 
TTNEW  = O.ODO 
TKI.EW  a U.ODo 
TENEW  a O.ODO 
TRHO  8 O.OLO 
TPRE  r O.ODO 
TVOL  a O.oDO 
TGAM  a O.ODO 
TENT  8 O.ODO 
TRHOC  a K.ol'O 
Pfl  50  I * 1»  >0 
xM ASS f I ) * O.ODO 
X I NEW ( I ) « O.ODO 
XKNFWCI)  a O.ODO 
XENFH(I)  « O.Ocn 
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r — 


63 

XRhP(I)  » 0.OD0 

6tt 

XPRE(I)  * O.Cf'O 

65 

XG AM ( I ) * 0.00c 

66 

XFNT(I)  * C.ODO 

67 

XRH8CCI)  * O.CDO 

68 

XRAPCT71)  ■ O.OPO 

69 

XVEL ( !♦ ' ) a O.oPP 

70 

XAPEAfl*) ) * O.OPO 

7l 

?0 

XVCU1)  * O.OOo 

72 

C 

73 

c 

the:  various  suhhaTions  are  perfprhed. 

76 

JSHL  * HATL'PCC?) 

75 

XPAPCl)  s PALH(I) 

76 

XVELCl)  ■ VNt«(D 

77 

XAREA(l)  * ArEa(1) 

78 

II  * l 

79 

UP  58  I * 2,  NP 

«0 

IHLD  a I 

Pi 

ISHL  * RATEPC(I) 

82 

IF  (ISHL  .LT.  1)  GO  TP  300 

83 

IF  (ISHL  .GT.  5)  GO  TP  301 

86 

IF  (ISHL  .E5.JSHL)  GO  TO  56 

85 

JSml  s IShl 

86 

XR AD ( I !♦ 1 ) s RAOM(I-I) 

87 

XVEL  ( 1 1 + 1 ) * Vf;EW(I-n 

88 

XAPEACII+l)  a APEA(I-t) 

89 

XENEW(II)  a XlNER(II)  ♦ XKKEW(II) 

90 

II  a II  ♦ 1 

•n 

56 

XKMEH(Il)  ■ XKNEW(II)  t KN£k ( I ) *LAM ( I ) 

92 

XIliFW(II)  a XlNfW(II)  + ENEW(I)*LAH(I) 

93 

XRhP(II)  a XPHO(H)  t PHEW(T)*LAM(I) 

9« 

XPRE(tl)  » XPReCII)  ♦ PNEW(I)*LAH(I) 

95 

XGAM(II)  a XGAH(II)  + GAMMaC (I ) *L AH ( I) 

96 

XENT(II)  * XeNT(II)  + ENTCCI)PLAM(I) 

97 

XRHPCCin  a XPhOCCII)  ♦ PhPC C 11 *L A M ( I ) 

98 

X VOL (II)  * XVOL(II)  ♦ LAH(I) 

99 

58 

XHASS(II)  a XHaSS(II)  ♦ PNEW(n*LAH(I) 

00 

XEKEW(H)  a XINEW(II)  ♦ XFlpw(II) 

»1 

Ic  « II 

02 

XPAD(IE*lI  a RaDNO.P) 

0 3 

XVEL(IE*1)  » VNEP(NP) 

no 

XAREA(IE+1)  = aREACHF) 

05 

L"  60  I = 1,  I I 

06 

TIMER  a TIME*  t XIMEwm 

07 

TKNFw  s TKf.FR  + XKMER(I) 

0 8 

TFMEr  s TEf.Ew  ♦ XEMfR  ( 1 1 

0 9 

TRMP  r TRHp  + XRHP(I) 

1 0 

TPRE  = TPPL  ♦ XPRE  (I) 

11 

TGAM  s T G A M + XGAH(I) 

12 

TFI.7  e TENT  ♦ XFMT(I) 

1 3 

TPhOC  S TFhOC  ♦ X R H P C ( I ) 

U 

TVOI  s tVPl  ♦ XV?L(n 

is 

6G 

TMASS  * TRASS  ♦ X«ASS(I) 

1 6 

TRhP  a TRHO/TVOL 

17 

TPRE  « TPPE/TVfiL 

1 8 

TGAH  a TGAM/TVCL 

19 

TFMT  a TENT/TVOL 

80 

TRrPC  a TRhO(./TvPL 

2 1 

PP  70  J * li  IE 

22 

XRHP(J)  = XPhOCJj/XVOLCJl 

23 

XPRE (J)  a XpPE (J) /X VPL CJ) 

^u 

XGAR(J)  a XGA“(J)/XVO( (J) 

i 
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XFNT(J)  = Xfcl.T  ( JJ/XVOl  ( J } 
XPhOCIJ)  = XFHCC(j)/Xv„LCj) 


' 


1?5 

l?b  7,t 


127 

C 

i?e 

C 

T^E 

LAYER-BY-LAYER  cross  TABULATIONS  ape 

PRINTEO. 

1 ?b 

II  * II  * 1 

1 Jr 

WRITE  ( b » 2d''9) 

131 

KRITEC6,  ?Ofb)  CXRADCn,  1 = 1,  II) 

1 32 

w R 1 7 E ( o » 2007) (XVELCI),  1 = 1.  II) 

1 J3 

WPIT[  (6»  ?00P.)  (XAREA(I) , 1 = 1,  II) 

l3u 

HP  I TP  (6#  1d()0) 

US 

WR I TE (b , lOOb)  TINCW,  (XINEw(J),  J=1 

* IE) 

13/, 

WP I TE  C b , 1007)  TKNf>,  (XKNEW(J),  Jsl 

» IP) 

t J 7 

WPITEfb,  1 0 OS)  TENE.W,  (XtNEWCJ),  Jsl 

. IF) 

1 38 

W P I T F.  f b » ur(>6)  TRhf,  (XRue(J),  J=l, 

IE) 

13o 

wPlTE(b»  0 on  7 ) TPPe,  f XPRf  { .1 ) » Jrl, 

IF) 

1 «o 

WRITE  (R,  oooR)  TGAf,  (XGAM(J),  j s 

1,  IE) 

t « 1 

WRITE  (6,  On l 0 ) TFNT,  (XfNT(J),  J * 

l,  It) 

10? 

W P I Tf  (b,  0 n 1 I ) TR^fC,  (XRhPC(J),  J 

* 1,  IE) 

103 

NP I Tt (b  * 0006)  TVOL,  (XVOLCJ),  Jsl, 

It) 

1 Oil 

WRITE  (Oi  0,.(I5)  TNASS,  (XNASS(J),  J = ) 

/ IP) 

1 os 

WRITE  (b,  0012) 

l"b 

RE  TURN 

107 

C 

1 os 

c 

POINT  A WARNING  IP  A*'V  SHFLU  PPIMTEP  IS 

LESS  THAI.  ONE. 

109 

3o  u 

WRITE  (b,  30o 1 ) I Si  L , IHLP 

uo 

Jo'i  1 

rflf'AT  ("'ISliL  = ',  u,  • USS  Than 

1 AT  CELL  10) 

151 

RETURN 

15? 

c 

I5J 

c 

print  a warning  ip  any  shell  pointer  is 

greater  than  FIVE. 

ISO 

301 

WRITE  ( b » 3011)  ISHL,  IHLD 

155 

1 

F«FMaT  {'CISHL  s IJ,  • GREATER  THAN  5 AT  CELL  10) 

15b 

RETURN 

157 

End 

u 
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SUBROUTINE  In  I T AL  (maTFp,  Gahma,  fr*C0N,  RaD,  VFL,  RHo,  pP£,  • ) 

I N I T A L SETS  UP  a p a T H E P GENERAL  HmltILaYER  FLUID  DYNAMICS  FPO- 

BL£m.  up  T 3 FIVE  Layers  OF  CONSTANT  (BUT  DIFFERENT)  DF'$Ity,  ppes- 
SUPF,  and  LTIFaPLY  VARYING  VELOCITY  CAt,  BE  INITIALIZED.  EACH  L A v t R 
can  have  A r if  F ERENT  EQUATION  OF  state  BUT  klTHJN  EACH  LAVEP  »Hp- 
COn,  GAMHa,  and  Trtc  MATERIAL  HATEP  APE  FIXED.  A VARIABLE  SPAClfG 
ACROSS  The  SHFLL  Is  ALLOWED  UNDEF  CONTROL  bf  the  pA°AmE  tEp  POX'S. 
PROVISION  IS  ALSO  made  for  a SINUSOIDAL  VELOCITY  PE°TURBAT  I 'f  . 


PARAMETER  NPT  s 2o? 


12 

REAL** 

RN, 

VN, 

PC, 

VO 

t J 

REAL*B 

mate rs, 

GAMMAS , 

PHBS, 

RhOCS 

1 u 

rfal*b 

R(  es, 

dELTar, 

HE  1.  T A V 

15 

RE  AL * A 

afc-i  , 

D R H (J , 

r-VEL » 

PC"S 

lb 

REAL** 

MATE  P (»'PT  ) , 

GAm,,A  (NPT ) , 

RH5C0N ( 

17 

REAL  ** 

pAD(NpT), 

VEL  (NPT  ) , 

Rho(npt)  , 

PRE (NPT ) 

1 A 

C 

1 0 

DATA 

n S H L L L « RN, 

VN  /l,  1. 

on-2o,  o.pro/ 

20 

DATA 

MCD( < DVEL, 

DRHO  / 1 - 

1. CD-2,  O.ODO/ 

2 1 

DATA 

lcells,  maters,  RHOCS 

/O,  1 . 0 1 D 0 , 0 

. o D o / 

NAMELIST  /SHLlNI/  NSHFLL,  RN,  VN,  NODE,  DP-r,  EVFL 
NAMELIST  /SHLOat/  LCEILS,  Rn,  VN,  ‘ aTFps,  GAMmas,  pm?S,  P«wS, 
PTES*  RhOCS 

format  ('cnamEliSt  /shldat/  data  for  layep1,  I?) 


initialize  the  loop  over  layeps. 

format  CPNAMILIST  /ShLINI/  DEFAULTS 
«R  I TE  (fc*  IOC  J ) 

•.RITE  (6*  SHLINI) 

REA'-CS*  SHLINI) 

format  f'CJAMFLIST  /ShLINI/  UPDATES  ...>: 
W P I T E ( 6 * t'U'J) 

WRITE  (b*  SHLlNI) 

RAt  ci)  - ph 
VEL(l)  1 VN 
1?  « 1 

L 9 100  ISHELL  = 1*  NSHELL 

pops  « l.CDO 

RO  s RAPCI2) 

VO  e V t L ( 1 2 ) 

REaD  IN  T«E  DATA  PER  The  SHELLS  one  At  a t^e. 
READ  (5  * SHLDAT) 

WRITE  (E,  ICO))  ISHELI 
aPITF  (6,  S H L r)  A T ) 

1 1 « T2  ♦ 1 

I?  * II  ♦ LCf  LLS  - 1 

SET  THE  SHELL  DATA  INTO  THE  OUTPUT  ARRAYS. 

(0  ?(iO  I s II*  I? 

PmCCOMCI)  * RHOCS 


GA'(Ma(I) 

MATER(I) 


Gammas 

maTlrs 


RMW(I)  s rk-os 

PPten  « PRES 

CELTAV  S DPLOATd  - I)  ♦ 1) 
rifLTAR  • Df.LT* v**PCWS 
pad(I)  * peltar*rn  ♦ r i .one 


1 ) /DPLOaT (12 


- deltar)»pp 


1 


63 

VEKI)  * 

DEL  T A V* VN  ♦ fl.OOO  - DELTAV)*VO 

tu 

200 

C out  ir.uE 

65 

C 

«>6 

t 00 

continue 

67 

N c 12  - 

1 

68 

►iP  e M ♦ 

1 

69 

C 

70 

c 

*00 

IN  A SINUSOIDAL  VELOCITY  PERTURBATION. 

7 1 

DO  3l)0  I 

s 2,  t.P 

72 

ARGI  s 3 

. 1015q26535eR79P0*DFLOAT(N0DE)*PADCI)/RAD(NP) 

73 

300 

VEL(I)  * 

VEL(I)  + DVEL*D5IN(ARGI) 

7 U 

FFTURM 

75 

END 
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Sl'e^PUT  INE  MAXMINU,  N,  AMAX,  IMAX,  AMIN,  JMJN) 


1 

2 c 

J C THE  maximum  (AMAX)  AND  MINIMUM  (AMIN)  OF  THE  N DOUBLE  PRECISION 

4 C VAUUCS  nr  VECTOR  A are  DETERMINED  and  the  CORRESPONDING  INDICES  IN 

5 C THE  APRAY  (AS  WELL  as  the  EXTREME  VALUES)  ARE  RETURNED.  THE  maxVaL 

6 C and  MINVAL  FUNCTIONS  COMPILE  AS  INLINE  VECTOR  OPERATIONS  BUT  ARE 

7 C RATHER  UNWIELDY  TO  USE  AND  DO  NOT  COMPILE  PROPERLY  UNDER  The  FX 

H C FORTRAN  COMPILER. 

9  C 

10  PARAMETER  NpT  s202 

11  REAL*?  A ( N) , AMAX,  AMIN 

12  INTEGER  IMAX,  i K l N 

1 3 C 

1 4 I«A X s MAXVAL  (A)  + 1 

15  I M I M s MINVALtA)  + 1 

Ifa  AMAX  s A(IHAX) 

17  AMIN  a A(IRJN) 

ia  return 

1 9 C 

2 0 C 

21  ENTRY  CIJBLOG  (DA,  DLNA,  N) 

22  C 

?3  C THE  CVECTORIZEO)  DOUBLE  PRECISION  LOGARITHMS  OF  THE  N INPUT  VALUES 

?4  C IN  REAL*8  ARRAY  (A  ape  COMPUTED  ANT  RETURNED  IN  THE  ARRAY  DLNA. 

25  C AS  "F  SUMMER  1°76,  THE  ASC,  FOR  SOME  OESCURE  REASON,  DID  NOT  have 

2b  C A DOUBLE  PRECISION  VECTORIZED  LOGARITHM.  THIS  ROUTINE  FILLS  THaT 

27  C PURPOSE. 

28  C 

29  RE  AL  *8  OA(N),  DLMA(N),  SCR(NRT),  DL (NPT ) 

30  REAL  SA(i.PT),  SLNA(NPT) 

31  C 

32  DO  10  I = 1,  '• 

33  $ A ( I ) = D A ( I ) 

34  10  SLNA(I)  s ALOGCSAd)) 

35  C 

ia  C*  20  I s 1,  N 

37  CLNA ( I ) = SLNA(I) 

38  2vj  SCRCI)  * DEXP  (DLi'A(I)) 

39  C 

99  r. " 30  I * 1,  N 

«1  SCRCI)  = OA ( I ) /SCR  Cl)  - l.ODO 

“2  L'L(I)  = SCR  ( I ) * ( 1 , ODC  - SCP  ( I ) * CO  . 5D0  - 0 . 1 666686600  *SC  R ( I)  ) ) 

93  30  OLNA(I)  s DLNA  (I)  + D L C ! ) 

99  RETURN 

95  END 


L; 
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Appendix  C 
STANDARD  TEST  #1 


1 AN  ALIABATIC  SOUND  WAVE  TEST 

2 

J &COI.TI  L 

4 M.NL 

5 &SHLINI 

o &END 

7 &SHLDAT 

6 tCLLLSslj,  RUbJ.OPO#  G AMHA Ss 1 . 4D0 , RHBSel.RDO,  PRESbI.ODO, 

9 FPWSsl.OOO, 

10  SEND 

1 1 
12 

13  THE  TIMINGS  F-  OR  ADimC,  EVEN  THOUGH  IT  IS  VECTORED,  ARE  A l8T 

u slower  than  cprpesfpnping  explicit  and  eulerian  calculations,  the 

lb  APINC  PACKAGE  IS  CONSTRUCTED  FOR  flexibility  AND  accuracy,  not 

lb  RAW  SPEED.  SEVERAL  CALCULATIONS  USING  STANDaPD  TEST  #1  HAVE  BEEN 

17  FFTFOPhED  T9  CONSTRUCT  TuE  FOLLOWING  TIMING  TABLE.  NCTEl  HOST  OF 

lo  THl  COMPETING  in  AC  INC  OCCURS  IN  THE  TRANSCENDENTAL  FUNCTION 

1 9 EVALUATIONS  F OP  Thl  EQUATION  OF  STATE. 

?o 


2 1 

22 

?3 

C 4 

UCELLS 

cpu  TJHE 

PC*  call 

cpu  time  per 

ITERATION 

1 0 

0.0082 

SEC 

0 . 0n23 

SEC 

?tl 

15 

C.009R 

StC 

0.0026 

SEC 

27 

?fc 

3o 

o . 0 1 3 2 

SEC 

0.0036 

SEC 

?9 

3o 

5o 

0.0175 

SEC 

0.0046 

SEC 

J 1 

It 

100 

0.0276 

SEC 

0.0073 

SEC 

33 

34 

2('o 

O.o«75 

SEC 

0.011° 

SEC 

V 
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Appendix  E 

STANDARD  TEST  #3 


1 A LINUS  SIMULATION 

2 

3 SCONTRL 

4 OTMINni.oD-7,  DTMAXpI .OD-4, 

5 EPSPobo.4500,  tPSV0«0.«5D0, 

b MAXSTPbIOOI.  IPRINTb50, 

7 ALPMA*2<  MCAMPelO, 

6 (END 

9 ISHLINI 

lb  NSHE  LL«5,  DVLLbO.ODO, 

1 1 SEND 

12  SSHLDAT 

1J  LCELLSbS,  PNbIO.OLO#  VNbo.OD-3,  maTEPS«i.01CO,  GAf'»'ASei  .666700, 

U RH0Sb1.2D-3,  PHOCSen.ODO,  PRE$el.oD6, 

15  SEND 

lb  SSHLCAT 

17  LCELLSb2,  RN«15.0D0,  VNao.oD-3,  MATEPSb2.0200,  GAMMASe2,oortCPO, 

IS  RHOSal .2D-3,  PhOCSbo.ODO,  PPES»1,006, 

19  SEND 

20  SSHLDAT 

21  LCELLSbIS,  PN«30.0Co,  VN«0.0£>-3,  haTEPS*3.03£>P,  GAVNAS*0, 5</0C'P</, 

22  PPWSb2.PD0,  RHOCS«l.oO<),  PRES«1.oD6,  PHOS«  1 . OOOOUCOOO 1 ODo , 

23  SEND 

24  SSHLDAT 

25  LCELLSbS,  PNbjs.oDO,  VNso.oD-3,  maTEPSb4 , 0^00 , GAHMAS»P.50«0Do, 

2b  PPKSbi .500,  RHOC5e7.8Do,  PHES*l .vDb,  RHPS*7.80000Pooo78ro, 

27  SEND 

2b  SSHLDAT 

29  LCELLSs5,  PNB40.0DO,  VNbo.OD-3,  maTEPSbS.oSD" , GAMHASei .4C00DC, 

3u  RHOSb  1 . 2D»  1 , PHBCSbo.PDO,  PPESsl.0D8,  POWSbI.ODO, 

3 1 SEND 
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LCELLS*  5,p"  * J5.tr'iP0O(>n()Bii«ac,V'i  * '!toooooo<,Onooon<,nc»MAtFPS*  «. <159999900900990, 
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T. 799999999999099, 
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